{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step-by-step Guide: Constructing a model with Phydrus\n",
    "\n",
    "*Authors: R.A. Collenteur & M. Vremec*\n",
    "\n",
    "This notebook is part of a manuscript that is currently being prepared (spring 2020): \n",
    "\n",
    "*R.A. Collenteur, G. Brunetti, M. Vremec & J. Simunek (in preparation) Phydrus: a Python implementation of the Hydrus-1D model.*\n",
    "\n",
    "---\n",
    "\n",
    "This Notebook presents the basics steps to create a Phydrus model simulating water flow through the vadose zone. In the presented example, the workflow is divided into several steps, to demonstrate the usage of Phydrus methods:\n",
    "1. Import the Phydrus package\n",
    "2. Create the basic model\n",
    "3. Add processes and materials\n",
    "4. Add soil profile\n",
    "5. Add Root Water Uptake\n",
    "6. Add atmospheric boundary condition\n",
    "7. Add observations nodes\n",
    "8. Simulate the model\n",
    "9. Post-process the results\n",
    "___\n",
    "\n",
    "### 1. Import the Pydrus package"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    " \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import phydrus as ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Create the basic model & add time information\n",
    "In the code block below a `Model` instance is created. The path to the Hydrus-1D executable has to be set at this stage. Phydrus will check the executable, and raise an Error if it is not found. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Folder where the Hydrus files are to be stored\n",
    "ws = \"output\"\n",
    "exe = os.path.join(os.getcwd(), \"../../hydrus\")\n",
    "\n",
    "# Create model\n",
    "ml = ps.Model(exe_name=exe, ws_name=ws, name=\"model\",\n",
    "              time_unit=\"days\", length_unit=\"cm\")\n",
    "\n",
    "ml.add_time_info(tinit=0, tmax=730);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Add processes and materials\n",
    "In this step a model for the water flow is selected and top and bottom boundary conditions, using the `ml.add_waterflow` method. After that, we can use the `get_empty_material_df` method to obtain an empty DataFrame to define our soil hydraulic parameters for the different soil materials. In this example, the model contains to soil materials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead tr th {\n",
       "        text-align: left;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th colspan=\"6\" halign=\"left\">water</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>thr</th>\n",
       "      <th>ths</th>\n",
       "      <th>Alfa</th>\n",
       "      <th>n</th>\n",
       "      <th>Ks</th>\n",
       "      <th>l</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3382</td>\n",
       "      <td>0.0111</td>\n",
       "      <td>1.4737</td>\n",
       "      <td>13</td>\n",
       "      <td>0.5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3579</td>\n",
       "      <td>0.0145</td>\n",
       "      <td>1.5234</td>\n",
       "      <td>50</td>\n",
       "      <td>0.5</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  water                                 \n",
       "    thr     ths    Alfa       n  Ks    l\n",
       "1   0.0  0.3382  0.0111  1.4737  13  0.5\n",
       "2   0.0  0.3579  0.0145  1.5234  50  0.5"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ml.add_waterflow(model=0, top_bc=3, bot_bc=4)\n",
    "\n",
    "m = ml.get_empty_material_df(n=2)\n",
    "m.loc[0:2] = [[0.0, 0.3382, 0.0111, 1.4737, 13, 0.5],\n",
    "              [0.0, 0.3579, 0.0145, 1.5234, 50, 0.5]]\n",
    "ml.add_material(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Add Profile information\n",
    "We develop the linear function of potential root water uptake distribution  $S^*_{p}(z)$ vs depth, following Hoffman and van Genuchten. \n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "S^*_{p} = \\left \\{\n",
    "\\begin{aligned}\n",
    "&1 && \\text{for} && z>L-r_1 \\\\\n",
    "&\\frac{z-[L-(r_1 + r_2)]}{r_2} && \\text{for} && L-r_1 \\geq z \\geq L-(r_1 + r_2) \\\\\n",
    "&0 && \\text{for} && L-r_1 z < L-(r_1 + r_2) \n",
    "\\end{aligned} \\right.\n",
    "\\end{equation} "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x11b03e210>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANAAAAGoCAYAAAA3omtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU9dn/8fdNQIMGkIAbBGUpSqVVREFcUDaLWpVfF1RqW5f2R12qPq2tj+CGrbVKbfu0tRtaCiiudUPFB4XWKlpEVFRwY4sCrkAChp1wP3/MSQgIycSZyZzvmc/runI5c9ZvYm7O/Tlzzom5OyLy+TTL9wBEQqYCEsmACkgkAyogkQyogEQyoAISyYAKKABmdoCZVZlZUfT+aTP7fhPt+2tmtjTa/+FmNt/MBkTzxpjZnU0xjrhSATUhMzvOzJ43s9VmtsrMnjOzPg2t5+7vuXuJu1ensY8xZrY5+oWvjPZ3dAbDvgX4YbT/V9y9p7s/ncH2EkUF1ETMrDXwGPAHoBToCFwPbMzB7u519xJgb2Am8KCZ2U7G1DyNbR0IzM/y+BJDBdR0DgJw97vdvdrd17v7k+7+GoCZNTOzq83sXTP72MwmmVmbaF5nM/M0f+FruftmYCKwH9DOzM6Njnq/NbOVwJhd7dfMdjezKqAIeNXMFkVjKTezITvbn5n1i454lWb2ak2rl2QqoKbzDlBtZhPN7GQza7vD/HOjr4FAV6AEuDWTHZrZ7tE2l7r7imjyUcBiYF/gF7var7tvjI5iAIe5e7cG9tUReBy4gdQR9ifAA2a2dybfQ9ypgJqIu68BjgMcuA34xMymmNm+0SJnA79x98XuXgWMAs5q7FEncoaZVQJLgSOAr9WZ9767/8Hdt7j7+izu99vAVHef6u5b3f0pYA5wyucYfzBUQE3I3d9093PdvQz4EtAB+J9odgfg3TqLvws0J3WkaKz73H0vd9/H3Qe5+0t15i3dYdls7fdAYHjUvlVGBXwcsH9jBx8SFVCeuPtbwARShQTwPqlfwhoHAFuAj7K96x3eZ2u/S4E7osKt+drT3W/KYKyxpwJqImbWw8wuN7Oy6H0nYAQwK1rkbuBHZtbFzEqAG0mdTduS46Fla793AqeZ2VAzKzKzYjMbUPP9JpUKqOl8SirAv2Bma0kVzjzg8mj+eOAO4BlgCbABuKQJxpWV/br7UmAYMBr4hNQR6ack/HfMdEOdyOeX6H8dRHItuAIys5PM7G0zW2hmV+Z7PFLYgmrhoosp3wFOBJYBLwIj3P2NvA5MClZoR6C+wMLoQ79NwD2kgqtIXnyeT7nzqSPbfxC4jNSZrVpmNhIYCVBcXHxEhw4dmm50khWfbEhd97p3cfa7o8WLF69w96xdXhRaATXI3ccB4wC6devmixY9m+cRSWOd+dcFANz7g+5Z37ZZx3cbXip9obVwy4FOdd6XRdNE8iK0AnoR6B59ar4bcBYwJc9jkgIWVAvn7lvM7IfANFL3qYx393pv9rrvvplNMjbJnk8+KQbC+H8X1GnsxtpZBtq8uZply6rYsKHBu6MlQ8XFRZSVldCiRVGj1stxBnrJ3Y/M1vaCOgJlw7JlVbRq1ZbOnduyk7ucJUvcnZUrK1i2rIIuXdrkezg5E1oGytiGDdW0a6fiyTUzo127tok/0if+CLRjH921a1cqKtbmaTSFZ+3aDY3OMspAMbGzDPTmm6v44hcPytOICs+bb77DF79Y2qh1QspABdfCxUFJScO/UN///gW88cabANx4483bzTvmmBM+1z7Ky8v50pcOT3OUjdO580GsWLGi4QUTpiBbuFWrqvI0mm0aGsPYsbfULnfjjTdzwQUX18577LHH0/oedlymsnId1dVbc/L9b926lYqKtTRrVrzd9KS3cIkvoDPOOG6792++uYrS0pJdLN10SktLePrpfzNmzA20b9+OefPmc8QRvbnzzgmYGQMGnMgtt9zEP/7xIOvXr2fQoIH07HkIkydPpKSklKqqVVRVVTFs2DeoqKhk8+bN3HDDGIYNO327fdS1Zs0egPPf//1Tnn/+P3Ts2IFHHnmAli1bsmjRIi6++DI++WQFe+yxB7fd9id69OjBo48+xg033MSmTZto166UyZMnsu+++7Jy5UpGjPgOy5e/z9FH98PMaNt2z8/s86OPij/z/6AhD0Qt3BlnZK3TqnXmmdndXuILqD7XP/o2b7z/aVa3eUiHVlx32sFpL//KK3OZP/8VOnTowLHHDuC5557nuOOOrZ1/002/4NZb/8zcuS9+Zt3i4mIeeuh+WrduzYoVK+jXrz+nn35avWcYFyxYyN1338Ftt/2ZM874Fg888BDf/va3GDnyIv7yl1vp3r07L7wwm4suuox//nMaxx13LLNmPYuZcfvt4xk79tf8+tdjuf76GzjuuGO59tqrePzxqfztb39v3A8qIQq6gOKgb98jKStLPXejV69DKS9/d7sCqo+7M3r0NTzzzEyaNWvG8uXv89FHH7Hffvvtcp0uXTrTq9dhABxxRG/Ky9+lqqqK55+fxfDh36pdbuPG1BOHly1bzplnns0HH3zIpk2b6NKlMwDPPDOTBx+8F4CvfvUU2rbd8TmRhSHxBVRfBrrk2I452We6+WTNmvU0a9a8dvnNm7dSWVnFqlVVbNlSzerV62rn7bjNVauquOuuu1i+/EOeemo6LVq0oFevXnz44Sp2261kp+tUVq6jefMWtdM3btzC2rXrWLFiDW3atOGf//zXZ/Zx4YWXcuGFF3LyySczc+ZMxo4dy6pVVVRXb6Wyctv43F0ZKIninIFat25JixZFteMpLm5BSUkxpaUlNG9eRJs2e1BaWkKLFi1o1Wp3WrRosd361dUbKSvbn333bcu//vU0S5cuZa+99qjd3s4yUFFRs9rpe+yxG1u3bqZz5w507dqFGTOmMXz4N3B3XnvtdQ477FDWrq2iR49ulJaW8NBDD9C8eWq8Awcez9Spj3L11aN44on/pbKysiAzkE5jB2DkyO9x6KFHcPbZ52w3/eyzRzBnzst8+cu9mTRpMj16pJ+9djR58gT+9re/c9hhR9KzZy8eeeRRAMaMuZrhw0dwxBH9aN++Xe3y1113Nc888yw9e/biwQcf4YADDvjc+w6ZPkiVnEr6B6mJb+Hi+jlQoUh6BtIRSHIq6UcgZSCRDKiFk5xSCxcwtXD5pxZORHZJBcRKYEUWv1bmdLRz577K1KlPNHq9999/n29+86x6l6nvdoeJE++ge/dD6N79ECZOvKPR+0+qgs9ApaXZbmE9pxlr5swXmDt3Lv369U97nS1btlBc3Jpx426vd2y7ut2hoqKC6677OTNmzMDMGDRoEP37D2SvvfZqcN9Jz0CJL6CGL+XZkPV91nepUHl5OSeddBr9+h3F88//hz59juS8877Lddf9nI8//pjJkyfSt28fZs9+kcsuu5wNGzbQsmVL/v73cXTp0oWxY8eyfv165sx5kVGjruDUU0/hkkt+xLx589m8eTNjxlzNsGGnM2HCJB588GGqqtZSXV3NxIm3c+qpX2PevFcoLy/nO985n7VrU7e233rr/3DMMUd/5lKfGtOmPc7QoUPo1i31TMuhQ4fwwgvPMWJEw9fFJP1SnsQXUBwtXLiI+++/m/Hjx9GnzzHcdde9zJz5L6ZMeZQbb7yZhx/+Bz16HMyzz/6T5s2bM336DEaPvpYHHriXn/3sWubMeYlbb/0dAKNHX8OgQQMYP34clZWV9O17LEOGDAbg5Zfn8tprcygtLaW8vLx2//vssw9PPTWV4uJiFixYwIgR32XOnP/scrzLly+nU6dtD4QtKytj+XI9EBYKoIAabuGyv8+G2qQDDzyQjh07U1m5ji98oTtHHXU0FRVr6dSpK4sWLWHVqiqWL/+AK6+8ksWLF2NmbN68mVWrqqiq2sCGDZtr9zF16jQeemgKN9/8awDWrVvPa6+9RVXVBo4//nhgN1atqtquPVuzZg1XXHEF8+bNo6ioiEWLFn1mmbrWrdvEhg0ba6evX78J92Zptapq4QIXtxZuzZo9aNmyuHaZli13p127NpSWlrBmTQnuWyktLeHHP/4VQ4cO5tJLH6K8vJwBA75CaWkJJSXFFBe3qF2/qKgZDz98HwcfvP2FpG+/PZ/S0ja1y9Vtz37/+99ywAEdueeeSWzdupXi4tbR/nfewnXv3oWnn/537fSVKz9mwIAT0rqqPektnM7CxdTq1avp2DF1v9KECdvOerVq1YpPP932L//QoSfyhz/8iZrP8155ZW5a295///1o1qwZd9wxmerq+p/dNnToiTz55HQqKiqoqKjgySenM3ToiZ/n20ocFRDZfsBidrZ3xRWXM2rU1Rx+eF+2bNn2F+cHDjyBN954k169+nDvvfdzzTWj2bx5M4ceegQ9e/bimmvGNLjtiy76ARMn3slhhx3JW2+9zZ577lnv8qWlpVxzzWj69DmGPn2O4dprr6I0F71vgBJ/JcIvf/nL7aZ17dqVrl275mlEhWfx4sUsXry4Uev85Z1UBrrgoOy312eeeWZWr0RIfAHpUp780qU8IrJLKiCRDCT+NLZuZ8ivpH8OpAwkOaUMJCK7VPAt3C9OOZSqVZ9kbX8lpXtz1dTXsra9Hb3++ut8+OGHnHhi4z7I/OCDDxg1ahQTJkzY5TLvvfceI0aM4LnnnvvMvOHDhzNnzhz69evH3XffnfZ+k97CJb6AGrqUJ5vFU7O9XD64ccmSBcyZ8xJnnvm1tNfZsmULpaXdmTLlH/Uut6tLeQBGj/4p69at469/vb1R358u5ZGsKi8vp0ePL3Puud/noIN6cvbZ5zB9+gyOPXYA3bsfwuzZqYfIz579IkcffTyHH96XY445gbfffptNmzZx7bU/4957/1F7JcLatWs5//yR9O17LIcf3pdHHpkCwIQJkzj99K8zaNBQBg8+abub5crLy+nffxC9ex9F796p2yoaMnjwIFq1apW7H0ygEn8EiqPQbmeQXUt8AeXjNHbSbmeosWbNejZvrm7Uz08ZKHD5eLh80m5nqLHjw/DToQwkeRGn2xlk1wq+gFq33yeW24vT7QwA/fsPYvjwbzFjxr8oK+vKtGlPZvLtJUbir0TQ7Qz5pdsZAqZLefJPl/KIyC4l/izczk5jr1z5ab1/yVqyw90Tfxq74Fq4JUtW06pVW9q1a6siyiF3Z+XKCj79tIIuXdo0at2QWrjEH4F2VFZWwrJlFXzyyYp8DyXxiouLKCvL/x90zqWCK6AWLYoa/S+iyK4kvoBC6KNle8pAMbGzDCTxF1IG0mlskQyogEQyoAwksaMMFBPKQGFSBhIpEGrhJHbUwsWEWrgwqYUTKRAqIJEMKANJ7CgDxYQyUJiUgUQKhFo4iR21cBkws18BpwGbgEXAee5eGc0bBXwPqAYudfdp9W1LLVyY1MJl5ingS+5+KPAOMArAzA4BzgJ6AicBfzKzoryNUoQYFpC7P+nuNU8SnAWURa+HAfe4+0Z3XwIsBPrmY4wiNeKegc4H7o1edyRVUDWWRdO2Y2YjgZEA7du3D6KPlu2FlIHyUkBmNh3YbyezrnL3R6JlrgK2AJMbs213HweMg1QGauyDzSX/Qnq4fF4KyN2H1DffzM4FTgUG+7azHMuBTnUWK4umieRN7DKQmZ0EXAGc7u7r6syaApxlZrubWRegOzA7H2MUqRHHDHQrsDvwVPTgw1nufoG7zzez+4A3SLV2F7t7g3+XI4Q+WrYXUgaK3edA2aTPgcKkz4FECkQcW7isCqENkO2phYsJtXBhUgsnUiBUQCIZUAaS2FEGiglloDApA4kUCBWQSAaUgSR2lIFiQhkoTMpAIgVCLZzEjlq4mFALFya1cCIFQgUkkgFlIIkdZaCYUAYKkzKQSIFQAYlkQBlIYkcZKCaUgcKkDCRSINTCSeyohYsJtXBhUgsnUiBUQCIZSHQGcvcg+mjZXkgZKNEFZGbMuPo7+R6GNFJF/58BMGPStXkeScPUwolkINFHIHdn8A135HsY0kgL3km1cIOHZv//3bgs/43HRBeQWrgwqYUTKRAqIJEMJLqFUwYKkzJQTCgDhUkZSKRAqIBEMpDoFk4ZKEzKQDGhDBQmZSCRApHoI5BauDCphYsJtXBhUgsnUiBUQCIZSHQLpwwUJmWgmFAGCpMykEiBUAGJZCDRLZwyUJiUgWJCGShMykAiBSLRRyC1cGFSCxcTauHCpBZOpECogEQykOgWThkoTMpAMaEMFCZlIJECkegjkFq4MKmFiwm1cGFSCydSIFRAIhlIdAunDBQmZaAsMLPLgVuAvd19hZkZ8DvgFGAdcK67v9zANpSBAqQMlCEz6wR8BXivzuSTge7R10jgz3kYmsh2YllAwG+BKwCvM20YMMlTZgF7mdn+eRmdSCR2LZyZDQOWu/urqa6tVkdgaZ33y6JpH+yw/khSRyjat2/H4D8qA4VGGagBZjYd2G8ns64CRpNq3z4Xdx8HjAPo1q2bKwOFJ6QMlJcCcvchO5tuZl8GugA1R58y4GUz6wssBzrVWbwsmiaSN7Fq4dz9dWCfmvdmVg4cGZ2FmwL80MzuAY4CVrv7BzvfUu32dBo7QGrhcmMqqVPYC0mdxj6voRV0GjtMauGyxN0713ntwMX5G43IZ8X1NLZIEGJ9BMqUMlCYlIFiQhkoTCFlILVwIhnY5RHIzL6exvob3H1qFscjEpT6WrjbgEcAq2eZ40mdXo4lZaAwJSUDPeHu59e3spndmdXRZJkyUJgSkYHc/dsNrZzOMiJJ1uBZODMrAr4KdK67vLv/JnfDyg61cGFKSgtX41FgA/A6sDWre88xtXBhCqmFS6eAytz90JyPRCRA6XwO9ISZfe77c0SSLJ0j0CzgITNrBmwmdVrb3b11TkeWBcpAYUpaBvoNcDTwenRFdDCUgcIUUgZKp4VbCswLrXhEmkI6R6DFwNNm9gSwsWZiCKexRXItnQJaEn3tFn0FQxkoTInKQO5+fVb32ISUgcKUqAxkZk+Z2V513rc1s2m5HZZIGNJp4fZ298qaN+5eYWb71LdCXKiFC1OiWjig2swOcPf3AMzsQLZ/5G5sqYULU0gtXDoFdBUw08z+TepD1P5Ej84VKXTpnET4XzPrDfSLJv2Xu6/I7bBEwlDfLd37ufuHAFHBPFbfMnGkDBSmpGSgqUDvBtZPZ5m8UQYKU1Iy0GFmtqae+QbUN18k8XZZQO5e1JQDyQW1cGFKSgsXPLVwYQqphdODFUUyoAISyUBaLVz0ZJ592f6pPO/teo14UAYKU6IykJldAlwHfMS2p/I4EPsHjSgDhSmkDJTOEegy4GB3X5nrwYiEJt1bulfneiAiIarvUp4fRy9rbul+nMBu6VYGClNSMlCr6L/vRV91b+nW7QySM4nIQDW3cpvZcHe/v+48Mxue64GJhCCdkwijgPvTmBY7auHClIgWzsxOBk4BOprZ7+vMag1syeoockQtXJgS0cIB7wNzgNOBl+pM/xT4US4HJRKK+jLQq8CrZnYXqVsXepA6efC2u29qovGJxFo6GehE4K/AIlKF1MXMfuDuT+R0ZFmgDBSmRGSgOn4DDHT3hQBm1g14HIh9ASkDhSmkDJTOlQif1hRPZDGpHCRS8NI5As0xs6nAfaQy0HDgRTP7OoC7P5jD8YnEWjoFVEzqSuwTovefAC2B00gVVGwLSBkoTInKQO5+Xlb32ISUgcKUqAxkZgeZ2Qwzmxe9P9TMrs790ETiL50W7jbgp6ROZePur0WfDd2Qy4Flg1q4MCWqhQP2cPfZZlZ3mi7lkZxJVAsHrIg++3EAM/sm8EFORyUSiHSOQBcD44AeZrac1J97/HZORyUSiHTOwi0GhpjZnkAzdw/mQ1RloDAlIgPVuaV7x+lAGLd0KwOFKaQMlM4t3QcDfYAp0fvTgNm5HJRIKNK5pfsZoHdN62ZmY0hdTCpS8NI5ibAvUPf+n03RtNhTBgpTIjJQHZOA2Wb2UPT+/wETsjqKHFEGClNSMhAA7v4LM3uC1B8XBjjP3V/J7bBEwpDWw+Xd/WXg5RyPJevUwoUpaS1csNTChSmkFk5/H0gkAyogkQzEsoWL/ibRxUA18Li7XxFNHwV8L5p+qbtPq287ykBhUgbKgJkNBIYBh7n7RjPbJ5p+CHAW0BPoAEw3s4PcvbqebSkDBUgZKDMXAje5+0YAd/84mj4MuMfdN7r7EmAh0DdPYxQBYngEAg4C+pvZL4ANwE/c/UWgIzCrznLLomnbMbORwEiA9u3bMfiPauFCoxauAWY2HdhvJ7OuIjWmUqAfqYtY7zOzrulu293Hkbp/iW7durlauPCE1MLlpYDcfciu5pnZhcCD7u6kLiHaCrQHlgOd6ixaFk0TyZs4ZqCHgYGQeiIQqb+Kt4LU7RRnmdnuZtYF6I5uq5A8i2MGGg+Mjx6jtQk4JzoazTez+4A3SD3U5OL6zsCBTmOHKqQMZKnfzWTq1q2bDylalu9hSCPNjjJQ32ezn4HGLdj0krsfma3txbGFEwmGCkgkA3HMQFmjDBSmkDJQogtIl/KEKaTPgdTCiWQg0UcgtXBhUgsXE2rhwqQWTqRAqIBEMpDoFk4ZKEzKQDGhDBQmZSCRAqECEslAols4ZaAwKQPFhDJQmJSBRApEoo9AauHCpBYuJtTChUktnEiBUAGJZCDRLZwyUJiUgWJCGShMykAiBUIFJJKBRLdwykBhUgaKCWWgMCkDiRSIRB+B1MKFSS1cTKiFC5NaOJECoQISyUCiWzhloDApA8WEMlCYlIFECkSij0Bq4cKkFi4m1MKFSS2cSIFQAYlkINEtnDJQmJSBYkIZKEzKQCIFQgUkkoFEt3DKQGFSBooJZaAwKQOJFIhEH4HUwoVJLVxMqIULk1o4kQKhAhLJQKJbOGWgMCkDxYQyUJiUgUQKhApIJAOJbuGUgcKkDBQTykBhUgYSKRCJPgKphQuTWriYUAsXJrVwIgVCBSSSgUS3cMpAYVIGiglloDApA4kUiEQfgdTChUktXAbMrBfwF6AY2AJc5O6zzcyA3wGnAOuAc9395Qa2pRYuQGrhMjMWuN7dewHXRu8BTga6R18jgT/nZ3gi28SxgBxoHb1uA7wfvR4GTPKUWcBeZrZ/PgYoUiN2LRzwX8A0M7uFVIEfE03vCCyts9yyaNoHdVc2s5GkjlC0b9+OwX9UBgqNMlADzGw6sN9OZl0FDAZ+5O4PmNkZwN+AIelu293HAeMAunXr5spA4QkpA+WlgNx9lwVhZpOAy6K39wO3R6+XA53qLFoWTRPJmzhmoPeBE6LXg4AF0espwHctpR+w2t0/2NkGRJpKHDPQ/wd+Z2bNgQ1EeQaYSuoU9kJSp7HPa2hD+hwoTCFlIHP3rG4wTrp16+ZDipblexjSSLOjDNT32exnoHELNr3k7kdma3txbOFEghHHFi5r1MKFKaQWLtEFpEt5whTSaWy1cCIZUAGJZCDRLZwyUJiUgWJCGShMykAiBUIFJJKBRLdwykBhUgaKCWWgMCkDiRSIRB+B1MKFSS1cTKiFC5NaOJECoQISyUCiWzhloDApA8WEMlCYlIFECoQKSCQDiW7hlIHCpAwUE8pAYVIGEikQiT4CqYULk1q4mFALFya1cCIFQgUkkoFEt3DKQGFSBooJZaAwKQOJFIhEH4HUwoVJLVxMqIULk1o4kQKhAhLJQKJbOGWgMCkDxYQyUJiUgUQKhApIJAOJbuGUgcKkDBQTykBhUgYSKRCJPgKphQuTWriYUAsXJrVwIgVCBSSSgUS3cMpAYVIGiglloDApA4kUCBWQSAYS3cIpA4VJGSgmlIHCpAwkUiASfQRSCxcmtXAxoRYuTGrhRAqECkgkA4lu4ZSBwqQMFBPKQGFSBhIpECogkQwkuoVTBgqTMlBMKAOFSRlIpEAk+gikFi5MauEaYGbDgTHAF4G+7j6nzrxRwPeAauBSd58WTT8J+B1QBNzu7jelsR+1cAFSC9ewecDXgWfqTjSzQ4CzgJ7AScCfzKzIzIqAPwInA4cAI6JlRfIqL0cgd38TUkeIHQwD7nH3jcASM1sI9I3mLXT3xdF690TLvtE0IxbZubhloI7ArDrvl0XTAJbuMP2onW3AzEYCI6O3G8eljnYhaQ+syPcgGiH7411wJQBzs7rRWgdnc2M5KyAzmw7st5NZV7n7I7nar7uPA8ZFY5jj7kfmal+5ENqYQxxvNreXswJy9yGfY7XlQKc678uiadQzXSRv4vY50BTgLDPb3cy6AN2B2cCLQHcz62Jmu5E60TAlj+MUAfJ3GvtrwB+AvYHHzWyuuw919/lmdh+pkwNbgIvdvTpa54fANFKnsce7+/w0djUuN99BToU25oIer7l7NrcnUlDi1sKJBEUFJJKBxBaQmZ1kZm+b2UIzuzLf49mRmf3KzN4ys9fM7CEz26vOvFHRuN82s6F1puf9ezKzy83Mzax99N7M7PfRmF4zs951lj3HzBZEX+fkYayXRD/j+WY2ts707P183T1xX6RONCwCugK7Aa8Ch+R7XDuM8StA8+j1zcDN0etDovHuDnSJvo+iOHxPpD5KmAa8C7SPpp0CPAEY0A94IZpeCiyO/ts2et22Ccc6EJgO7B693ycXP9+kHoH6El364+6bgJpLf2LD3Z909y3R21mkPtuCOpczufsSoOZypjh8T78FrgDqnnkaBkzylFnAXma2PzAUeMrdV7l7BfAUqesbm8qFwE2euiwMd/+4zniz9vNNagF15LOX/nTcxbJxcD6pf8Vh12PP6/dkZsOA5e7+6g6zYjle4CCgv5m9YGb/NrM+0fSsjjdu18IlSjqXM5nZVaQ+85rclGPbmfrGC4wm1XbGRgPjbU6qfewH9AHuM7Ou2R5DUguovkuCmow3cDmTmZ0LnAoM9qhBJ4+XM+1qvGb2ZVJ54dXoCvoy4GUz61vPeJcDA3aY/nRTjDca84XAg9HPdbaZbSV14Wt2f775DNI5DJDNSYXWLmwLhD3zPa4dxngSqSsu9t5hek+2D7mLSQXc2HxPQDnbTiJ8le1PIsyOppcCS0idQGgbvS5twjFeAPwsen0QqR0MAXUAAAGtSURBVPbMsv3zTeQRyN23fM5Lf5rSraT+Jz4V/as+y90v8OxfzpRrU0mdiVsIrAPOA3D3VWb2c1LXMULql3lVE45rPDDezOYBm4BzPFVNWf356lIekQwk9SycSJNQAYlkQAUkkgEVkEgGVEAiGVABiWRABRQAM+scfZ7R0HITzGyJmV2QhX2eGV3W/1im20oyFVDy/NTd/5LpRtz9XuD7WRhPoqmAwlFkZrdFN4c9aWYtG1rBzPaNbtZ7Nfo6JjqavRUdrd4xs8lmNsTMnotufOvb0HZlGxVQOLoDf3T3nkAl8I001vk98G93PwzoDdRcmvIF4NdAj+jrW8BxwE9IXXUtaVIBhWOJu9c87fYloHMa6wwC/gzg7tXuvrrOtl53962kimpGdJ3Y62luVyIqoHBsrPO6msxuRam7ra113m/NcLsFRwWUbDNI3dpM9Gdi2uR5PImjAkq2y4CBZvY6qbZPf1Mpy3S4DoC7lwNfqvP+ljTX+4idPxij7rbO3dV+pGE6AiXLauDn2fogFfgTUJHxqBJMN9SJZEBHIJEMqIBEMqACEsmACkgkA/8HKUl6ePnNiEIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 216x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define loop for potential root water uptake distribution proposed by Hoffman and Van Genuchten\n",
    "def z_loop(z, r1 = 10, r2 = 20):\n",
    "    if z > -r1:\n",
    "        return 1\n",
    "    elif z < -(r1 + r2):\n",
    "        return 0\n",
    "    else:\n",
    "        return(z+(r1+r2))/r2\n",
    "\n",
    "bottom = [-30, -100]  # Depth of the soil column\n",
    "nodes = 150  # Dictretize into n nodes\n",
    "ihead = -500  # Determine initial pressure head\n",
    "\n",
    "profile = ps.create_profile(bot=bottom, dx=1, h=ihead, mat=[1,2])\n",
    "profile[\"Beta\"] = profile.apply(lambda row: z_loop(row[\"x\"]), axis=1)\n",
    "ml.add_profile(profile)\n",
    "ml.plots.profile(figsize=(3,6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Add atmosphere boundary conditions\n",
    "Atmospheric boundary condition can be added easily by reading in a CSV file using Pandas `read_csv` method and adding it to the model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "atm = pd.read_csv(\"../data/atmosphere.csv\", index_col=0)\n",
    "ml.add_atmospheric_bc(atm, hcrits=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Add root water uptake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.add_root_uptake(model=0, p2h=-1500, p2l=-1500, poptm=[-25, -25])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7. Add observation nodes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Number of the node -- > write script to get node closest to desired depth\n",
    "ml.add_obs_nodes([-30, -60])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 8. Write hydrus input files & run hydrus \n",
    "Before we can simulate, we write all model information to files. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Successfully wrote output/SELECTOR.IN\n",
      "Successfully wrote output/PROFILE.DAT\n",
      "Successfully wrote output/ATMOSPH.IN\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "CompletedProcess(args=['/Users/Raoul/Projects/phydrus/examples/phydrus_paper/Ex1/../../hydrus', 'output', '-1'], returncode=0)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ml.write_input(verbose=False)\n",
    "ml.simulate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9a. Plot results\n",
    "Plot pressure for soil column at the end of the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([<matplotlib.axes._subplots.AxesSubplot object at 0x11b94d410>,\n",
       "       <matplotlib.axes._subplots.AxesSubplot object at 0x11b9b0a50>],\n",
       "      dtype=object)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAC2CAYAAADUZoqGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXwU9f348dd7jxzkIJCEJCQc4ZQjQRAUUbEKgtiKtkWrqKDWavlW1Ep/KtXWo+3ja4utX7WoLVVbLZZapbS1qPTyAKEIcslluEmEAAFysjk/vz9mExbIsQnZnT3ez8djH5mdmZ15z2R23jvX+yPGGJRSSkUvh90BKKWUspcmAqWUinKaCJRSKsppIlBKqSiniUAppaKcJgKllIpyLrsDaK+0tDTTt29fu8NQEWrt2rVHjDHpdsxbt20VSK1t22GXCPr27cuaNWvsDkNFKBHZa9e8ddtWgdTatq2nhpRSKsppIlBKqSiniUAppaJc2F0jaE5tbS2FhYV4PB67Qwk5cXFx5OTk4Ha77Q4lNOxZDutfhyufhLhku6NpW+UReP9Ju6MIT84YcPj81nW4weECEXC6ISbResUmQZfukNADknuGx3bRySIiERQWFpKUlETfvn0REbvDCRnGGEpKSigsLCQ3N9fucELD/tWwfiFc+mB4fOGry+Gzt+yOIgwZqK8F0+B9a6ChFhrq2v5oYgZk5kOvC+CcL0PG0MCGGgIiIhF4PB5NAs0QEVJTUzl8+LDdoYSO2hPW36697I3DX91z4cHddkcRWYw3SdRWWom2uhyqSqDiEJTuh8Ofw4H18J+fwH9+bCWES74HgybZHXnAREQiADQJtEDXy2lqq8CdcOopAxVdRMAVY73iu7U8XvlB2PxnWPUCvH4dDP4yTH0OElKDF2uQ6Lehk4gIc+bMaXr/1FNP8dhjj7VrGomJie0af+3ateTl5TFgwADuuecetG0JP9RWgTve7ihUOEjKhLGzYPZamPg47PgnLLgMDm+3O7JOp4mgk8TGxrJ48WKOHDkStHnOmjWLBQsWUFBQQEFBAe+++27Q5h22aj2aCFT7ON1w8X1w21Lr1OJrX7OOFiKIJoJO4nK5uPPOO3n66afPGLZnzx4uv/xy8vPzmTBhAvv27QNg9+7dXHjhheTl5fHII4+c8pl58+YxZswY8vPzefTRR8+Y5oEDBygrK2Ps2LGICDNmzGDJkiWBWbhIUucBV6zdUahwlDMabn4LThyDP9xoXWeIEBFzjaDR43/bzJYvyjp1mkN7JvPo1cPaHO873/kO+fn5PPDAA6f0nz17NjNnzmTmzJm8/PLL3HPPPSxZsoR7772XWbNmMWPGDObPn980/rJlyygoKGD16tUYY5g6dSoffvgh48ePbxqnqKiInJycpvc5OTkUFRV1wtJGtoZaDyca3LjrGohxhf7vIE9tPZuKSu0OI2w5RHA7BcG6ViYCsS6HdQepw0Gsy0FinIuEGBdOhx/X07Ly4asvwBszYMUzMP57AV6C4Ii4RGCn5ORkZsyYwbPPPkt8/MnTDytXrmTx4sUA3HLLLU2JYsWKFbz11ltN/R988EHASgTLli1j5MiRAFRUVFBQUHBKIlAds6e4hNJjtfSoqCY7JfRPER0s9XDdiyvtDiMqJMW66J4YQ0ZSHDnd4slNSyC/Vwpj+najS4zPrnLoNTD0Wvjgp1Z32kD7gu4kEZcI/PnlHkj33Xcfo0aN4rbbbvNr/Obu6jHGMHfuXO66664WP5ednU1hYWHT+8LCQrKzs9sfcJRxN9RQg4vUhBi7Q/FLRnIcv//mBXaHEbbqGhqoqz95E0W9MdTUNWCAhgaDp7aeiuo6yj11lJ6opaSyhuJSD6t2lbB4nXWEHeNyMH5gGtPOy2HS0EwcDoGr5lkXj99/Eqa9ZNPSdZ6ISwR26969O9dffz0vvfQSt99+OwDjxo1j0aJF3HLLLSxcuJBLLrkEgIsuuohFixZx8803s3DhwqZpTJ48mR/84AfcdNNNJCYmUlRUhNvtpkePHk3jZGVlkZyczKpVq7jgggt49dVXmT17dnAXNhzVV1PniCXO7bQ7Er/Exzi5eGCa3WFEpTJPLev3Hec/2w+xdNMB/rn1EMN6JvPja4czsncPGPNN+Pg5+NJcSBtgd7hnJfRPkoahOXPmnHL30HPPPccrr7xCfn4+r732Gs888wwAzzzzDPPnzycvL++U8/uTJk1i+vTpTReSp02bRnl5+Rnzef7557njjjsYMGAA/fv3Z8qUKYFfuDDnqK+hwREeRwPKXslxbsYPSufRq4fx8UMTeOaGczlWWcN1L67kNx/twlx4NzhjYfkv7A71rEm43Xs+evRoc3rN9q1btzJkyBCbIgp9un5OOviT4XxOb8Y/vLTZ4SKy1hgzOshhAc1v2yq0lJ6o5cE3N/Lu5oPcN3Eg93lesGpXzdkO8Sl2h9eq1rZtPSJQUcXRUEODU28fVR3TNd7N8zeN4rrzcvi/fxaw1D3RuiV582K7QzsrmghUVHGZWowmAnUWHA7hya/nc9ngdO77UPB0GwTrFrb9wRCmiUBFFbepscoTK3UWnA7h59efS2piLL+tugiK1sCRHXaH1WGaCFRUcZlacMXZHYaKAN0TYnh86jBeLbWe92H73+0N6CxoIlBRJYZaLTGhOs0VQzPo238w2+hL3dbmb0AIB5oIVPSor8NJA6JHBKqTiAjfv2oI79WNwlG42mpRLgwFNBGIyJUisl1EdojIQ80M/7aIbBKR9SKyXETCtikgO8pQP/zww/Tq1avdn4tWDd5GaRxuTQSq8wzP7kpx5uU4aKBu+3t2h9MhAUsEIuIE5gNTgKHAjc3s6F83xuQZY84FfgaE7ZMZdpShvvrqq1m9enXQ5hfuPJ4qABwxempIda4Jl19BsUnh4Nq/2R1KhwTyiOB8YIcxZpcxpgZYBFzjO4IxxrdMaAIQXk+3+Qh2GWqAsWPHkpWV1fkLE6E8J6xE4NT2CFQnu+ycDDa4RpB4YJXVFGaYCWStoWxgv8/7QuCM6lki8h3gfiAGuPys5/rOQ3Bw01lP5hSZeTDlyTZHC2YZatV+1Scajwj01JDqXA6HQO4lpOz4gMO7N5Leb4TdIbWL7ReLjTHzjTH9gQeBR5obR0TuFJE1IrImlBti9y1D7WvlypVMnz4dsMpNL1++HLDKUN94441N/Rv5lqEeNWoU27Zto6CgIEhLEblOeCoBcMV0sTmSk8Jl21ZtGzLuKwBsXxV+t5EG8oigCOjl8z7H268li4AXmhtgjPk18Guw6rG0Olc/frkHUrDKUKv2q/FYF4vdsaFzRNCubVuFtF79hlDs6AG7P7I7lHYL5BHBJ8BAEckVkRjgBuCvviOIiG+LDl8Gwv5nr28Z6kaNZaiBZstQN/ZvNHnyZF5++WUqKioAqzWyQ4cOBWsRItbJRKDXCFRglGaMZWjNRvaXVNodSrsELBEYY+qAu4H3gK3AG8aYzSLyhIhM9Y52t4hsFpH1WNcJZgYqnmAKVhnqBx54gJycHKqqqsjJyWn37arRpsZ7aig2PsHmSFSk6j5kPN2lgtVrw+tuPi1DHQV0/VhWLX2NsavvpvC6peQMu6jZcbQMtTorh7bB8xfwq27f4657f2B3NKfQMtRKAXU11l1DcXpEoAIlbRAeZyJdj6yj3FNrdzR+00SgokZdtTcRdNFEoALE4cDTYyQjpIDlBeFTbkITgYoaddXWxeIu8VqSQwVO0sBxDJL9rN6+z+5Q/BYxiSDcrnUEi66XkxpqvLWGYvSuIRU4zt7n4xRD2Y5Vdofit4hIBHFxcZSUlOhO7zTGGEpKSoiLC5375u3UWHQOLTGhAin7PAAyyjdzoPSEzcH4J5APlAVNTk4OhYWF6JOZZ4qLiyMnJ8fuMEJD7QnqceDUFspUIMV3oya5N8OP7ebjHSV8/bzQ//5FRCJwu93k5ubaHYYKcVJbRbXE0qWZp7mV6kzunFGMKPuYn+88EhaJICJODSnlD0fdCWpET5OpwJOeI8nmENt27bE7FL9oIlBRw1F3ghqHJgIVBD3PBSC1bGtYXCfQRKCihrPeQ71TE4EKgiyrDHWe7GbNnmM2B9M2TQQqargaPNQ79Y4hFQTx3TDdchnh2sPavZoIlAoJxhhiGjw0aMP1KkgkawTnuvbxyZ6jdofSJk0EKiqcqK2nCx4a3FpeQgVJZh4Z9QcoPHAw5OsOaSJQUaHsRB3xVIMmAhUs3usEg9nLxsJSm4NpnSYCFRVKT9SSIB4csZoIVJBk5gEw1LGXdftC+zqBJgIVFY5X1dAFD864JLtDUdEiMQMS0hnbpYh1+47bHU2rNBGoqGAlgmrcWnlUBYsIZOYz3LmPdfuPh3QtNE0EKiqUV1bikgZi4vWIQAVRZh5Z1bspr6xi/9HQfbBME4GKCifKrUPz2IRkmyNRUSUzD6epY4AUsW5/6F4n0ESgooKnwvoSxiak2ByJiiqZ+QCMcO8P6esEmghUVPBUWF9CR3xXmyNRUSW1P7jiuTjpABsKNREoZauaSu+XMFZPDakgcjghYyjDnfvZ/EUZNXUNdkfULE0EKirUVXkTQZwmAhVkGcPp6Smgpq6e7QfL7Y6mWZoIVFQwnjKrQ48IVLBl5hFTU0omR1kfoqeH/GqhTES+DAwDmip2GWOeCFRQSnU2qfYmgji9RqCCzPuE8dguRWzYf5xbxvaxOaAztXlEICIvAt8AZgMCXAeE3pIo1YKaugbctRXWm1h9jkAFWcYwAMYnH2T9/tA8IvDn1NA4Y8wM4Jgx5nHgQmBQYMNSqvOUVFbTVSqpdXYBp9vucFS0iU2CbrnkOfez83BFSFYi9ScRND4OVyUiPYFaICtwISnVuQ6WekiRcupiu9kdiopWmXlkV+/AGNhUFHqVSP1JBG+LSAowD/gU2AP8wZ+Ji8iVIrJdRHaIyEPNDL9fRLaIyEYR+ZeI6Ckn1emKyzx0owK6aCJQNsnMp0vFXhI4wYb94ZkIfmaMOW6MeQvr2sA5wI/b+pCIOIH5wBRgKHCjiAw9bbR1wGhjTD7wJvCz9gSvlD8OlnroJhU4E1LtDkVFq8zhAFyacpgNIXidwJ9EsLKxwxhTbYwp9e3XivOBHcaYXcaYGmARcI3vCMaY/xhjqrxvVwE5/oWtlP+Ky6tJkQrciWl2h6KilffOoS8lHwzJJ4xbvH1URDKBbCBeREZi3TEEkAx08WPa2cB+n/eFwAWtjP9N4B0/pqtUuxSXeuguFUiX7naHoqJVcjbEdyPfvY8DpR6KyzxkJIdO+9mtPUcwGbgV61f6L3z6lwHf78wgRORmYDRwaQvD7wTuBOjdu3dnzlpFgeLSSpKohBBMBLptRwkRyBhOr4odAKzff5zJwzJtDuqkFk8NGWN+Z4y5DLjVGHOZz+saY8xiP6ZdBPTyeZ/j7XcKEZkIPAxMNcZUtxDLr40xo40xo9PT0/2YtVInlR07ggMD8aGXCHTbjiKZ+XQ5/jmxjoaQu07gzzWCFSLykoi8AyAiQ0Xkm3587hNgoIjkikgMcAPwV98RvKecfoWVBA61M3al2lRb30Bd2UHrTWIPe4NR0S1rBFLnYWL68ZB7sMyfRPAK8B7Q0/v+c+C+tj5kjKkD7vZ+divwhjFms4g8ISJTvaPNAxKBP4nIehH5awuTU6pDDhz3kGq8DYIkhc6huIpCWVbbBJd1PcjGwlLqG0Kn6Up/ag2lGWPeEJG5YO3gRaTen4kbY5YCS0/r90Of7ontCVap9tp7tJIeeBNBYoa9wajoljoQXPGMcO6lonooOw9XMCgjNEqe+HNEUCkiqYABEJGxQOg9EaFUM/YdrSJdvJurJgJlJ6cLMoeTU10AwLp9odN0pT+J4H6sc/v9RWQF8CpWATqlQt6+o1VkOo5jYhIhNtHucFS0yxpBXMkWusY5Quo6QZunhowxn4rIpcBgrGcJthtjQq9qklLN2HW4krGxFYheKFahIGsE8slvmJh5IqTaMParPQKsp4T7escfJSIYY14NWFRKdZLtB8vJcZdCol4oViEgawQAE7p+weK9sZR7akmKs78irj/tEbwGPAVcDIzxvkYHOC6lzlpFdR37jlaR0XAIumr1EhUC0oeAM4Y8x26MIWQK0PlzRDAaGGqMCZ17nZTyw+fF5bioI6nmEHTra3c4SoErBjLzyKrYgsjlrN17jIsH2l8Dy5+LxZ8Belytws62A+X0lBLENEA3rXCuQkTPUbiKN3JOehc+DZE7h1orOvc3rFtGk4AtIrIaaCoBYYyZ2tJnlQoF2w+WMSimxHqToolAhYjsUfDJAib1LeOVz2toaDA4HNL25wKotVNDT2HdJfRT4Fqf/o39lAppm78oY0JyGZSjRwQqdPQcBcBF8ft5xpPLjhB4sKzFRGCM+QBARNyN3Y1EJD7QgSl1Njy19WwsLOW7vUuh0mWVAVYqFKQNhJhEzmkoAHJZvfuo7YmgxWsEIjJLRDYBg71NSTa+dgMbgxeiUu23Yf9xauobGOA8BF17gcNpd0hKWRxO6DmSpJKN9EiKZfXuo3ZH1OqpodexGor5X8C3veFyY4z9kSvVisYvV1rVDuhxegupStksZzTy8XNc1C+BlbuPYoxBxL7rBK21R1BqjNljjLnRGLPX56VJQIW81XuOkpcRi/PYTsgYZnc4Sp0q53xoqGNy9wMcLPNQeOyEreH4c/uoUmGlrr6BT/ceY0pGKZgGyNAjAhVicsYAMMphtVhm9+khTQQq4ny67ziVNfWMSyq2evTQIwIVYhLToVsu6cc3kNLFzapdJbaGo4lARZylmw4Q43Iw1FkIrjjo3s/ukJQ6U6/zkcJPGNu3Oys1ESjVeRoaDO9tPsilg9KJOfwZpA+26sArFWpyxkBFMZOyrWsE+0qqbAtFE4GKKOsLj3Og1MOXh6VC4RrropxSoajPOAAucn8OwMc7j9gWiiYCFVHe2XQAt1OY2O0Q1FY2fdmUCjnpQyC+Gz2OriU9KZaPd9p3ekgTgYoYdfUN/H3jAS4ekEbigf9aPftcZG9QSrXE4YDe45B9HzOufyof7zxCg00N2msiUBFj2ZZivij18I0xvWHvx5A6AJK0nWIVwvqMg6O7uKKX4UhFDVsPltkShiYCFTF+89Eu+qR24Ypz0mDfx3paSIW+vtYR6yXu7QC8v/2wLWFoIlARYe3eY3y67zi3X5SL8+B68JRC30vsDkup1mXkQWwyXYtXMTQrmQ8+10SgVIe9tHwXyXEupp2XA58tBmcMDJxkd1hKtc7pgr4Xw673uXRwOp/uPUaZpzboYWgiUGHv8+Jy3v3sINMv6EOC2wGb/wz9J0B8it2hKdW2fl+CY3uYnHWCugbDioLg30aqiUCFtYYGw8N/3kRyvJs7x/eDwtVQ/gUM/5rdoSnln36XATC8Zh3JcS7+ufVQ0EPQRKDC2pufFvLJnmN8f8oQuifEWKeFXHEweIrdoSnln7SBkNQT1+4PuPycHvx7WzF19Q1BDUETgQpbRytr+N+lWxnTt5t1baDWA5sXw8ArINbeFp+U8psI9L8Mdn/AFUPSOVZVy9q9wW3UXouwqLD1k79vpdxTx4+vzbMa/167ECoPw5hv2R1ap6mtraWwsBCPx2N3KCElLi6OnJwc3G633aF0joFXwPqFXNZlN26n8I8txVzQLzVosw9oIhCRK4FnACfwG2PMk6cNHw/8H5AP3GCMeTOQ8ajI8ftVe3nr00K+c1l/BmcmQX0drPg/yB4NuePtDq/TFBYWkpSURN++fW1twSqUGGMoKSmhsLCQ3Nxcu8PpHP0ngMNFlz3/YFz/r/Du5oM8/OUhQfufB+zUkIg4gfnAFGAocKOInN5CyD7gVqxmMZXyywefH+bRv27mssHpfHfiIKvnZ2/B8X1wyRzrUDtCeDweUlNTNQn4EBFSU1Mj6ygpLtkqh7L9Xb6Sn0XhsRNsKCwN2uwDeY3gfGCHMWaXMaYGWARc4zuCtynMjUBwr4yosLX9YDl3L/yUgT0SeW76KFxOBzQ0wPKnrSJeg660O8ROp0ngTBG5TgZPgSPbubLnCWKcDv624YugzTqQiSAb2O/zvtDbr91E5E4RWSMiaw4ftufJO2W/Q2Uebv/tJ8THOHn51jEkxnrPbH7yGzi8FcZ/zyrkFUbCYdsWEebMmdP0/qmnnuKxxx5r1zQSExPbNf7atWvJy8tjwIAB3HPPPRhjTzG2oPLe6Za05z3GD0rj7xsPBK0IXVh8a4wxvzbGjDbGjE5PT7c7HGWDTYWlXDN/BUcra3hp5hh6psRbA47uhn8+CgMmwvCv2xtkB4TDth0bG8vixYs5ciR4DzrNmjWLBQsWUFBQQEFBAe+++27Q5m2bbn0hawRsXsLVI3pysMzDqt3BKU0dyERQBPTyeZ/j7adUu/xlfRHTXvwYhwhvzrqQvJyu1oCGBvjL3eBwwdXPRtS1gVDicrm48847efrpp88YtmfPHi6//HLy8/OZMGEC+/btA2D37t1ceOGF5OXl8cgjj5zymXnz5jFmzBjy8/N59NFHz5jmgQMHKCsrY+zYsYgIM2bMYMmSJYFZuFAz9FooWsPk7DoSY128tTY4u8xA3jX0CTBQRHKxEsANwPQAzk9FmPoGw1PLtvPC+zs5v293nr95FGmJsSdH+OQ3sHc5TP0ldO3QWcew8vjfNrPli84tUzy0ZzKPXj2szfG+853vkJ+fzwMPPHBK/9mzZzNz5kxmzpzJyy+/zD333MOSJUu49957mTVrFjNmzGD+/PlN4y9btoyCggJWr16NMYapU6fy4YcfMn78yTu9ioqKyMnJaXqfk5NDUVGU/IYceg3863HiCv7GV/LH89cNX/DENcNIiA3snf4BOyIwxtQBdwPvAVuBN4wxm0XkCRGZCiAiY0SkELgO+JWIbA5UPCq8fHH8BN/83Se88P5Opl/Qm9/fccGpSWD3R7DsEeuU0Mib7Qs0SiQnJzNjxgyeffbZU/qvXLmS6dOt33e33HILy5cvB2DFihXceOONTf0bLVu2jGXLljFy5EhGjRrFtm3bKCgoCNJShIHU/pCZD5+9xbTzcqiqqWfppgMBn21A04wxZimw9LR+P/Tp/gTrlJFSANTUNfDyit08+68CGozhR9cO55axfU4dqWgt/OEG6J4LX1sQNaeE/PnlHkj33Xcfo0aN4rbbbvNr/Obu7DHGMHfuXO66664WP5ednU1hYWHT+8LCQrKzI/+Ir8mIG+C973Nel2L6pSfw+up9XDe6V9ufOwthcbFYRYeVO0u46tmPePKdbVw0II1/fPfSM5NA8Rb4/dehSyrcsgS6dLcn2CjUvXt3rr/+el566aWmfuPGjWPRokUALFy4kEsusdqAuOiii07p32jy5Mm8/PLLVFRUANZpoEOHTi2ylpWVRXJyMqtWrcIYw6uvvso115xy53lky7seHC5kwx+4+YI+rNt3nM+KAvtMgSYCZbvCY1Xcu2gdNy5YRXVdPS/fOpoFM0bTq3uXU0c8vB1e+yo4Y2HGXyA5y56Ao9icOXNOuXvoueee45VXXiE/P5/XXnuNZ555BoBnnnmG+fPnk5eXd8r5/UmTJjF9+vSmC8nTpk2jvLz8jPk8//zz3HHHHQwYMID+/fszZUoUFRFMTLfa0tjwR74+MpN4t5PXVu4N6Cwl3O7PHT16tFmzZo3dYaizZIxh+Y4jvLpyL//aWozL4eDbl/bjfy4bQJzbefrIsP51WPo9cHeBW9+GHkMCEpeIrDXGjA7IxNvQ3La9detWhgwJzLKGu4heN1vfhj/eBDe8ztwtvVn8aSEfP3Q5qb7XydqptW1bi86poCrz1PLW2kJeW7WXXYcr6Z4Qw7cv7c9NY/uQ3fhsgK/qcnj7ftj0htX05NcW6JGAinyDroTkbFi9gG9e+XsWfbKP3328h/snDQ7I7DQRqIBraDBsKirljTX7+fO6Iqpq6jm3Vwq/uH4EV+VlnXkE0OiLdfDm7XBsD1z2sFVHyNHCuEpFEqcLzrsN/vNjBjgOcsWQDH63ci93Xdo/ILeSaiJQAeGprefjnUf4x5ZD/HtbMcVl1cS4HEwd0ZMZF/YhP6eVZiRLC+H9J2H9QkjKglv/Dn3GBS94pULBeTPhg5/Cf1/krksfZtmWYhb+dy93ju/f6bPSRKA6zeHyav69rZh/bj3E8oIjnKitJyHGyfhB6UwYksGEc3rQLSGm5QlUHYWPfg6rFwAGLvg2jP9/emeQik6JPSD/G7DuNc679EHGD0rnxQ92Mf2CPifrbHUSTQSqw45W1rBh/3HW7TvGRzuOsH7/cYyBnl3juG50DhOGZDC2X3diXW2czvGUwepfwYpnoaYCRtwIX3oIUnoHZ0GUClUX32cdGa96nvuvuJdr56/gleW7mT1hYKfORhOB8ktNXQNbDpSxft8x1u8/zrr9x9lbUgWAQyAvuyvfnTiIiUMyGJKV1HaZ4Po62PU+bPgDbHsb6jxwzlfg8kcCdkeQUmEnbSAMnQqrF3DuuNlMGprBix/s5Bvn96JHUlynzUYTgTpDXX0D+4+dYFNRKeu8O/7NX5RRU2c1G9EjKZaRvVO4YUxvRvZOIS+7q/8XsIo3W7eCbvoTVBRDXIpVImLkzdBzZACXSnWUiHD//ffz85//HLDKUFdUVLSrFHViYmLTQ2T+ePjhh3n11Vc5duxYuz4Xkb40F7b8FZb/grlXzWXS0x8w793tzLtuRKfNQhNBFDteVcPOw5XsPFzBrsOV7Dpcwa4jlewtqaS23nq+JM7tID87hVvH9eXcXimc2yuFrK5x/jcMUl9n3f2z+wPYsgQObrKqhQ6cDOfeaD044+r4vdEq8BrLUM+dO5e0tLSgzPPqq6/m7rvvZuDAzj0FEpZ6DLHKTqxeQO4Fs7j9olx+9eEubji/F+f16ZzrZ5oIIlxtfQP7jlaxq2mH793pH6nkaGVN03hup9AnNYF+aQlMHJJBv/QEhmYlMzgzCbezHQ+gNzRA8Wew+0PrtfdjqPE+OdpzFEyZZ7UbkBC8hrnV2fEtQ/2Tn/zklGF79uzh9ttv58iRI6Snp/PKK6/Qu3dvdu/ezfTp06moqDijPL1t9iQAAA3kSURBVMS8efN44403qK6u5qtf/SqPP/74GfMcO3ZsQJcp7HxpLny2GP7xA2Zf/Wv+tuEL5i7exNuzLyHGdfYFIjQRhLHqunqKS6s5WObhQOkJiss8HCj1nPxb6qG4vJp6n1aO0hJj6ZeewORhGfRLS6RfegL90hPp1S3eavaxvepq4PA22P9f61f/nuVw4pg1LHUA5F9nNSbf9xJICM6vyYj1zkPWEVVnysyDKU+2OVowy1CrZnTrAxd/Fz54ksRRM3jimuHc8eoafvmfHdx/xaCznrwmghBkjKGiuq5ph964Uz9Q5v3r3dmX+Pyib5QQ4ySzaxxZXeMZNyCNrK5x5KZZO/vctAS6xrs7HljVUevX/sFNJ1+Ht0NDrTW8ay8Y/GVrx597CST37Pi8VEjxLUMdH3/yCfCVK1eyePFiwCo33ZgoVqxYwVtvvdXU/8EHHwROLUMNUFFRQUFBgSYCf1x8H2z8I7z9XSZ+ewVfG5nNL/9dwPiBaYzue3aniDQRBEF1XT3Hq2opqajhWFUNRytP/m18We9rOVpZzbHKWmrqG86YTmpCDBnJcWR1jePc3ilkJceR2TXOu+OPIyM5jqS4s9jRN6r1QOl+OLTFu8P37vzLTpYGJjHT+jU58Arrb89RVlN7UVIS2hZ+/HIPpGCVoVYtcMfD1Gfhd1fDv3/E49c8wZq9x7h30Xr+Nvtiurf2jE4bNBG0Q01dAxXVdVR46iivrqXCU8fxE7Und+aVNRyt8v5t6q6lorquxWmmdHHTvUsM3RJiyE6JJy87me4JsXRPcHt3+vFkJsfRIzm25VIM7VVXY+3Uj+2F4/u8L2/3sb1QcfDkuOKEtEHWk72Zw62dfkaeVSFRRRXfMtS33347cLIM9S233NJsGeqbb775jDLUP/jBD7jppptITEykqKgIt9tNjx49bFmmsJM7HsbcAateIGnARH45/TymvbCSexet47e3nY/T0bEfYhGfCIwxnKit9+68rZ14RXUd5d6/FR5rR11eXUdlc8Obdvx1TbdPtqRLjJNuXWLonmDt2PulJ3rfu+mWEENqQswpw1Pi3R07L9+a+jqoOgIVh7yvgyd39o07/vIvwPgsiziha471ANfAiZDSx+pOHwzpQ8Ddefcrq/A2Z84cfvnLXza9f+6557jtttuYN29e08VisMpQT58+nZ/+9KenXCyeNGkSW7du5cILLwSs20p///vfn5EIHnjgAV5//XWqqqrIycnhjjvuaNftqhHtiiesmzAWf4v8uz7iiWuG8dDiTfzo7S08evVQ/+/o8xExZah/9u42PvuijApPLZXV9d6dubWTb/BjEWNcDpJiXSTGuUiMtV5Jjd1xLhJiXdbwWBeJce6mcVK6uOmeYO3cO+0X++ka6qHyCFQesu69rzjs7fa+fLurSoDTF1isSoYpva2LTim9T+7su/WBpJ5WkSulZajDSFSvmyMF8OvLILUf3PYOP162l98s383cKedw16XN1yKKijLUh8urKT1RS1Ksix5JcU07dN+d+cn3bhJinSTFur07eWfbZRDOljHW07Oe0mZex8/sd+LYyR1+Vcmpv+AbueKteiSJPaBbLvS64OT7hB4nu5NzwNXx84dKqRCTNhCmvWQ12frmN/n+da9SXF5N0fETGGPafVQQMYmgM5+yO4UxUFcNtVVQU+n9WwE1VVa3p6z5HXlzr8a7a1rijLGetI3rar269YFeY07dqft2xyTqxVmlotWgyTDlZ7D0ezj+/C2enrYAp8vdoVNDEZMIqK+D2kprZ11T1Uy3dwdeW3Vmd7PjVp7c+Zt6/2JwxZ3ciccmQ3w3606axn6nvFLO7Kfn4pVS7XH+t6wfqssexlXngWmvQEyXtj93mshJBL+9ynqoyV+ueGuFxSSAO8HqdneB+JyT3TEJ3uHNdDf+9d3x645cBUBHDvUjXbhd2wyocXdbZVqW/j+recubF7f7TEHkJIIxd8CQqd6duHen3Wy3dyevLV2pMBAXF0dJSQmpqamaDLyMMZSUlBAXpz+8mpz/LUhI7/Dp4shJBPnX2x2BUp0uJyeHwsJCDh8+bHcoISUuLo6cnBy7wwgtw67t8EcjJxEoFYHcbje5ubl2h6EiXCc/zaSUUircaCJQSqkop4lAKaWiXNiVmBCRw8DeFganAUeCGE5rQiWWUIkDQieW1uLoY4yxpaJeG9t2R4TK+gaNpTnBjqPFbTvsEkFrRGSNXXViThcqsYRKHBA6sYRKHIEWSsupsYRuHKCnhpRSKuppIlBKqSgXaYng13YH4CNUYgmVOCB0YgmVOAItlJZTYzlTqMQRWdcIlFJKtV+kHREopZRqp5BOBCJynYhsFpEGERl92rC5IrJDRLaLyGSf/ld6++0QkYd8+ueKyH+9/f8oIjHe/rHe9zu8w/v6EdcfRWS997VHRNZ7+/cVkRM+w170+cx5IrLJO59nxVtBTES6i8g/RKTA+7dbO9fRYyJS5DPPqzp7HfkZxzwR2SYiG0XkzyKSYtc6aSPOZpc9ErS2LZw2XsDXQUvbQzPj7fFuA+tF5MymBzs+/1aXsSPf+w7G0UtE/iMiW7z7snubGedLIlLq83/7YSBiaZUxJmRfwBBgMPA+MNqn/1BgAxAL5AI7Aaf3tRPoB8R4xxnq/cwbwA3e7heBWd7u/wFe9HbfAPyxnTH+HPiht7sv8FkL460GxgICvANM8fb/GfCQt/sh4KftnP9jwPea6d9p68jPOCYBLm/3TxuXw4510kqMLS57JLxa2hbsWActbQ/NjLcHSOvkebe5jGf7vW9HLFnAKG93EvB5M7F8CXjbzm0npI8IjDFbjTHbmxl0DbDIGFNtjNkN7ADO9752GGN2GWNqgEXANd5fmpcDb3o//zvgWp9p/c7b/SYwofGXaVu8410P/KGN8bKAZGPMKmP9519tYf6+cZ2tzlxHbTLGLDPG1HnfrgJaLQ1p0zppdtk7adrhIijroL3bQyfzZxk7/L1vD2PMAWPMp97ucmArkN3Z8zlbIZ0IWpEN7Pd5X+jt11L/VOC4z4bZ2P+UaXmHl3rH98clQLExpsCnX66IrBORD0TkEp95FDYTF0CGMeaAt/sgkOHnvH3d7T0Ef9nnNEpnrqP2uh3rF34jO9ZJc1pa9kjS3Lbgy451cPr24MsAy0RkrYjc2Unz82cZz+Z73yHe008jgeZa0LpQRDaIyDsiMiyQcTTH9jLUIvJPILOZQQ8bY/4S7Hga+RnXjZx6NHAA6G2MKRGR84Al7fmnGmOMiJxxG1drsQAvAD/C+kL9COtU1e3+zrM9/FknIvIwUAcs9A4LyDqJVqGyLbQVSyvbw+kuNsYUiUgP4B8iss0Y82FgIraPiCQCbwH3GWPKThv8KVb5hwrvdZ0lwMBgxmd7IjDGTOzAx4qAXj7vc7z9aKF/CZAiIi5v9vcdv3FahSLiAroCJW3F5R33a8B5PstSDVR7u9eKyE5gkHcevofGvvMvFpEsY8wB7+mSQ6fPy991JCILgLdPW67m5tnedeRXHCJyK/AVYIL3dE/A1kkHtbZOwkIHtwVfnbYOOrI9NDONIu/fQyLyZ6zTOmebCPxZxma/92c532aJiBsrCSw0xiw+fbhvYjDGLBWR50UkzRgTtDpE4Xpq6K/ADd4r/7lY2XM18AkwUKy7X2KwLgL91bsR/geY5v38TOAvPtOa6e2eBvy7pY32NBOBbcaYptMbIpIuIk5vdz9vXLu8pznKRGSs9zzkjBbm7xuXX7w7ykZfBT7zmW5nrSN/4rgSeACYaoyp8ukf9HXSimaXvZOmbbtWtgVfQVkHLW0Pp42TICJJjd1YF5ibi7m9/FnGjn7v28W7bb8EbDXG/KKFcTIbr0+IyPlY++WAJKUW2Xmluq0X1sZciPWLshh4z2fYw1h3BmzHe7eJt/9VWFfmd2Idojb274e1I9wB/AmI9faP877f4R3ez8/Yfgt8+7R+Xwc2A+uxDveu9hk2Gmsj3wn8kpMP86UC/wIKgH8C3du5jl4DNgEbsTburM5eR37GsQPrnOt67+tFu9ZJG3E2u+yR8GppWwB6AkuDuQ5a2R6aYvFubxu8r82dGUtzywg8gZWYOvy970AcF2Odqtvosy6uAr7duP8A7vYu/wasC+vjgr3t6JPFSikV5cL11JBSSqlOoolAKaWinCYCpZSKcpoIlFIqymkiUEqpKKeJIEKJVfWzM+7JVkpFOE0ESinbiUhFJ09vj4ikdeY0I5kmgsjmFJEFYtVBXyYi8XYHpJQKPZoIIttAYL4xZhhwHOspX6VClljmichnYjVY8w1vf4e3Bs82sRorWioi09qY3GwR+dQ7nXOCEH7Y0kQQ2XYbY9Z7u9diNRKjVCj7GnAuMAKrntc8bw2lr2Ftv0OBW4AL/ZjWEWPMKKyqrN8LSLQRQhNBZKv26a4nBKrNKtWGi4E/GGPqjTHFwAfAGG//PxljGowxB7EKJLalsdKn/ghqgyYCpVSkavwhpD+C2qCJQCkVSj4CviEiThFJB8ZjVQddAXzde60gA6udX9VJNEtGKGPMHmC4z/un7ItGKb/9Gev8/was8s0PGGMOishbwARgC1Z560+xmpdUnUDLUCulwoKIJBqrOcdUrKOEi7zXC9RZ0iMCpVS4eFtEUoAY4EeaBDqPHhEopcKWt53j3NN6P2iMec+OeMKVJgKllIpyeteQUkpFOU0ESikV5TQRKKVUlNNEoJRSUU4TgVJKRbn/D+sosEyRPwbiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x180 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ml.plots.soil_properties(figsize=(6, 2.5))\n",
    "#plt.savefig(\"../../figures/soil_properties.eps\", bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9b. Plot the drainage over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAACsCAYAAAAuVDhiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2dZ3hcxdWA37NFWnXJkmzLVW6AuzHGxhjTmykmJIAhEEooISEEQkJxgFDyEUgIEAgkQEiCQwkQOg6YYqqp7r132ZJtybYkS1pJK53vx70rrWSVlbRN0rzPc5+9Ze7M2SvNnntmzpwjqorBYDAYDLGGI9oCGAwGg8HQFEZBGQwGgyEmMQrKYDAYDDGJUVAGg8FgiEmMgjIYDAZDTOKKtgChJisrS3Nzc6MthqEbsXDhwkJVzY62HKHC9CFDpGmuD3U5BZWbm8uCBQuiLYahGyEiW6MtQygxfcgQaZrrQ1Ed4hOR00VkrYhsEJHbmrgeLyIv29e/FZHcyEtpMBgMhmgQNQUlIk7gCWAaMAK4SERGNCp2JbBPVYcCjwB/iKyUBkPXoLWXQYMhFommBTUR2KCqm1S1CngJOKdRmXOAWfb+q8BJIiIRlNFg6PQE+TJoMMQc0ZyD6gtsDzjOAyY1V0ZVfSJSDGQChRGRMAJ4q2vYW1bF/vJq9ldUUemrpcpXS3VNwGeNUuWrJSXexWmjepOW4I622B2ioqqGrXvL2HugigOVPg5U+lCFQ3unMKpvWrTF64rUvQwCiIj/ZXBVWyrZVeLliU82MLRnMkOzkxnSM5meKfF0tXfG6upq8vLy8Hq90Raly+HxeOjXrx9ud3C/YV3CSUJErgGuARgwYECUpbGo9NWwYkcJq3YWs7qglJ37K4hzOshOiaemVtm2t5xNe8ooKGlbJ3j4w3V89KvjSI6P3T+dqlJa6aOg2MvmwjK2FJaxpajM3i9v8TtfPXUQvzljeJf70Ysyrb4MBtOHtu8t5/VFOzhQ6as7l+JxMSQ7mSHZyYzsk8rxh2YzODs51PJHlLy8PFJSUsjNzTX/hyFEVSkqKiIvL49BgwYFdU80f+V2AP0DjvvZ55oqkyciLiANKGpckao+DTwNMGHChCaj327ac4C73l6JiCCACAjgEGH6uD6cM65vh78QQHF5Nf/3v1XMWVFAqd2R0xLc9O+RQJWvloVb9yEi9O+RwNFDMxmUmURWSjwZiW5SE9wkuJ3EuRzEOR24nQ7iXPan08HcNbu46ZWlLN62j6nDYserucpXyxuL83h/5S427TlAQYkXb3VtgzI9kuLIzUys+84Ds5LomRJPcryLJFvZ/mPeJv7+xWYG9EjkR5Nz2bm/gn/O28z8LXs59pBsfnXqodH4et2CYPrQhNweLL/7VHaXVrJh9wE27D7Axj3W5xfr9/DaojzunQ25mYmcNLwX5x3Rj+E5qRH9HqHA6/Ua5RQGRITMzEz27NkT9D3RVFDzgWEiMghLEV0I/LBRmbeBy4CvgfOAj7Wd4ddrVSn1+lAAVdT6YE9pJZ+u20PPFA+Th2S297sA1hvmZf/6ju17y/neuL6cNLwnY/qlk5PmCck/+zFDswDYUlgWEwqqtlZ5dVEej360nh37KxiclcTIvmmcPLwXPVPj6ZXqITczidyspKCGJe+dPor8/V7uenslryzIY1V+CQB90xP4y8cb+P74fgzKSgr31+qKBPMyGBQiQq9UD71SPUyx/x/9bN9bzidrdzN39W6e+3or/5i3mcMHpHPRxAGcPaYPCXHO9n+DCGOUU3ho63ONmoKy55R+DrwPOIF/qupKEbkXWKCqbwP/AJ4TkQ3AXiwl1i6G9kzhzeumHHS+1FvN9574kuteXMT1Jw4lxePmrDE5eNxt60y7Srz88JlvKKnw8cJVRzFxUI/2itosmcnxOAT2HKgKed1tZcPuUma+vpz5W/Yxrn869507iuMOye5Qx3Y4hMcuOpwH31/L2oJSrp46mEuOGkBZZQ2n/flzFm/bZxRU+wjmZbDD9O+RyKWTc7l0ci77yqp4ffEOXvx2K7e8uoz/m72Kq6cO5opjBsX08LQhtojqf4qqvgu82+jcbwP2vcD54ZQhxePmqR9N4LJ/fsc971hzxq8tzOP5qybhdAT3Y6uq/OS5hew9UMWLVx/F2P7pYZHV6RB6JMWzp7QyLPUHy5qCEs5/8mtcDuGP543h/CP6heyNMynexd3TRzY4V1OreNwOVu4s4fvjQ9JMt6K5l8FwtpmRFMeVxwzix1Nymb9lH09/vpGHPlzHrK+38p+rJzGsV0o4mzd0EUwsPmBoz2Q+v+UEvvvNSdwzfSRfbyrixpeX8Pt3VzPz9WV8uGpXi/d/s2kvS7bv5zdnDg+bcvKTnRJdBZW3r5zL/vkdSXEuZv9iKhdM6B/24RCnQxiUlcymPQfC2k5XRlXfVdVDVHWIqt4XqXZFhImDevDMZUfyxs+Oxldby72z2+Q8aAgh5513Hps2bWqxjNPpZNy4cYwdO5bx48fz1VdftVh+y5YtvPjii3XHy5cv5/LLLw+FuEZB+XE6hJ6pHi6dPJBLjhrAO0t38uxXW5i9NJ+fPr+QD1YWUFvb9PTXP+ZtpkdSHD8Y3y/scmanxLPnQHQU1K4SL5c88y3lVTXM+vFE+qYnRKztIdlJbNxTBsDO/RWUeKsj1rYhNBw+IIOrpw7mi/WFbNhdGm1xuh0rV66kpqaGwYMHt1guISGBJUuWsHTpUu6//35mzpzZYvnGCmr06NHk5eWxbdu2DstsBoMbISL83/dG85szhpPgdlLi9XHmY19wzXML+cmxg5l5xvAG5bfvLWfuml1cd/zQNs9btYfs5Hg27o6sJVFW6eOj1bv445y17C+v4t9XTuLQ3pEdohmSncz/luezubCM6Y/P48zROTzwgzHtrs//nU44rCepns69rqwzMePI/vz5o3U8/822g4ZyY5F73lnJqp0lIa1zRJ9U7jq75e9eVlbGBRdcQF5eHjU1Ndx5553ceuutLFiwgKysLBYsWMCvf/1rPv30U+6++242b97Mpk2b2LZtG4888gjffPMN7733Hn379uWdd97B7XbzwgsvcM45ViyEJ598ko0bN/Lggw8C8Oyzz7JgwQIef/zxBnKUlJSQkZEBWFMZt9xyC++99x4iwh133MGMGTO47bbbWL16NePGjeOyyy7jl7/8JWeffTYvvfQSt9xyS4eelbGgmiExzoWIkJbg5qVrjgJg4dZ9B5X778I8AC6c2P+ga+HAP8TXTmfGNrN9bzknPfQZN7y0BJdTePHqozhiYEZE2g5kcHYSqnD+k19R6vWxOr/9PxrlVT5mPP01N7y0hPtmrw6hlIbWyEqO54zROby2MI/yKl/rN3RT5syZQ58+fVi6dCkrVqzg9NNPb7H8xo0b+fjjj3n77be55JJLOOGEE1i+fDkJCQn873//A+DLL7/kiCOOAOAHP/gBb7zxRt39L7/8MhdeaPmgVVRUMG7cOA477DCuuuoq7rzzTgBef/31Osvqo48+4uabbyY/P58HHniAqVOnsmTJEn75y18CMGHCBL744osOPwdjQQVBv4xEzj28L99t3tvgfE2t8t8F25k6LJt+GYkRkSU7JZ6qmlqKK6pJT4wLa1ve6hp+8txCyqp8vHDVJCYPzsQRpONIqBnZx4owUXigir7pCWzaU4aqtmv+6y8fb2DFjhKykuN5Y8kObpt2GBlJ4X2Whnp+dNRA3lqykxe/3cZVU1seboo2rVk64WL06NH86le/4tZbb+Wss85i6tSpLZafNm0abreb0aNHU1NTU6fQRo8ezZYtWwDIz88nO9tanpKdnc3gwYP55ptvGDZsGGvWrGHKFMvL2T/EB/D1119z6aWXsmLFCubNm8dFF12E0+mkV69eHHfcccyfP5/U1IPXuvXs2ZOdO3d2+DkYCypI+qR72FXipSZgHmr5jmLyi738YHxoFvkGQ+9UDwDb91aEva2X529nVX4Jj1wwjilDs6KmnMByZLn77BHcedYIrpiSS2mlj5KKtr+BV9fUMuurLUwf24fnrpxIla+Wd5Z1vCMZgueIgRkcd0g2j360nrx95dEWJyY55JBDWLRoEaNHj+aOO+7g3nvvxeVyUVtrLYBvHIYpPj4eAIfDgdvtrntxczgc+HxWP0lISGhw34UXXsgrr7zCa6+9xrnnntvky97kyZMpLCxs0+Jav3wJCR2fozYKKkh6pnjw1Sr7yuvXIK2xh5nGhdlzL5Cx/S1LYuYby/DV1LZSuv3U1CqzvtrC2P7pnDyiV9jaaQuXTxnElccMopetpHeVtj1W2rK8YsqrajhjdG+G56TSNz2Bbzftbf1GQ8gQEe6ZPhIErpq1gKIoOf3EMjt37iQxMZFLLrmEm2++mUWLFpGbm8vChQsBeO2119pc5/Dhw9mwYUPd8bnnnstbb73Ff/7zn7rhvcasWbOGmpoaMjMzmTp1Ki+//DI1NTXs2bOHzz//nIkTJ5KSkkJpaUOnl3Xr1jFq1Kg2y9gYo6CCxL+4sCwgDtnq/BKS4pz0j9DwHljDjXFOByt2lLB2V/g8oWYv28mmwjJ+cmzsDcH4FVRBcUMFpaqtzk19s8mKlDVxkBU15PAB6SzZvh+AbUXlFFcY78BIkJuVxF8vHs/mwjLOeeJL5q7eFbF51c7A8uXLmThxIuPGjeOee+7hjjvu4K677uKGG25gwoQJOJ1td8g688wz+fTTT+uOMzIyGD58OFu3bmXixIl15/1zUOPGjWPGjBnMmjULp9PJueeey5gxYxg7diwnnngif/zjH+nduzdjxozB6XQyduxYHnnkEQA++eQTzjzzzA4/B1S1S21HHHGEhoP3lufrwFtn64od++vOnf/kV3ruE/PC0l5LLNm2TwfeOlvfX5F/0LWXvtuqv31zeYfq99XU6ol/+kRPffgzramp7VBd4WBrYZkOvHW2vvTd1gbn31m6QwfeOlvnNPFc/FzyzDd62iOf1R0/88UmHXjrbL3+xUU68NbZeuWz37VZHqzIJ1H/3w/VFq4+1BSLt+3TE/70iQ68dbae9shn+thH6/TrjYV6wFsdMRkas2rVqqi1HU7Ky8t10qRJ6vP5wtqO1+vVSZMmaXV103/Dpp5vc33IOEkESb0FVQNYin1Nfglnje0TcVkyk60J/f3lDd/2a2uVW19bDsBNpxxKWmL73KdnL9vJxj1lPPHD8VGdd2qOvhkJpHpczF29mxlH1kfeXrLNsoTWFpRy2sjeB91XXVPLgi37mHFkvcfl2WNy+OOcNby91JqH+sYM90WUcf3Tef/GY3ltYR4vzd/OQx+uq7uWlRxH/x6J9ElLICs5jqzkeLJT4slKjicrJZ7MpDjSEt0kx7li8v801khISOCee+5hx44dYc36sG3bNh544AFcro6rF6OggiQp3jKp/UN8O4u9lHh9UYnW7Pfe21/RMCbf5qKyuv28/eWkJbY9t5Kq8tdPNnJIr2SmjTr4Rz4WcDqEaaNyeHdFfoPzRWVVddebYlleMRXVNUwKiJPYM9XD2z8/hlJvNV9uKOKRj9bhra6JyJo2PyISp1bSzm6J2+ngwokDuHDiAPaXV7Fw6z7WFJSyfW852/aWs7qghD2llZR6m3aKcYgVsiw1wUVKvJsUj4vkeBfJAZ8p8f59N8nxLsb0S6NPCwvNtZ0eorHOaaedFvY2hg0bxrBhw5q8pm0cxjUKKkj8FpQ/F47fQWJ4hBesAiTFOXE55CALym9BABS1M6Ds1xuLWLurlAfPGxPTb6UDMhMp9foaKBO/R1hFVU2T9yzeZq1jOyK34Tou/6LjrUXW/fnF3rAFpRWRT4HLVXWLfTwR+DswNiwNdjLSE+M4aXgvThp+sGOOt7qGorIqCksr2VNaSVGZpbRKKqopqfusptTro6DEy4E9Pg54fZRW+qjyNXQoSnA7+d33RnHeEQdHf/F4PBQVFZGZmdkllVS0UDsflMfjCfoeo6CCJKmRk8SaAstBIdIRFcDygkpPdLO/0YS+f7IfYG9ZvYJaubOY619czFlj+3D50bn0aGHNz3/mbycj0c3ZURi6bAv+9B0lFdV1Cmpzoa2gqptRUNv30zc9gZ4pTXcQ/xv1zv0V4Yyafj8wR0Qew0okOA24IlyNdSU8bid90xPaFWKr0ldDWWUNB7w+9pZX8Yf31vDr/y6luKKaK49pmDyvX79+5OXltdm12tA6/oy6wWIUVJAkNbKgVueX0L9HAilRCpOTluBmdX4Jy/L2M6af5ea+ZPt+RvVNZcWOEgoDXHef/nwTmwrLeGzuej5YWUB2SjyH9krhjrNGNKizoqqGj1bt4tzxfSM6xNUe/Apqf0U1PVM9lHir675zcwpqybb9LS4J6JNuKa784vCl+lbV90XkWuBDoBA4XFULwtagAYB4l5N4l5MeSXEMyEzk31dO5PoXF3Pf/1Yxsk8qRw2uzwXndruDzvhqCC/GzTxIPG7rUVXaQwWr80s4rHf0soVmJsezeNt+pj/+JWCZz+t3lzJpUCZup9TNx+wu9fLu8nyumJLL9ScOZU1BKV+sL+QfX24+qM6P1+ymorqGs8bkRPS7tAe/gvK7hW8prJ9/8zYxxLentJId+ytaVFB+9/X8/eFbBC0idwJ/AY4F7gY+FZEQ+OMa2oLb6eBPF4wlNzOJX/xnMfvKuu0UYExjFFSQuB3Wo6quqUVV2b43rMNArRI4zFFthz7yVtfSJz2BHklxdYsfX/puO9U1yqWTc7nhpGH89eLxXDxpAKrWsEcgs5ftJDslnkmDOpZZOBKk2x6KxfY83OYABdWUBeVfHzWqb/OOIx63k6zkOLaHN7pBJjBRVb9W1aeA04Abw9mgoWmS4108dtHh7C6t5IVvt0ZbHEMTGAUVJA6H4BDw1SglFT6qamrpmRIfNXl+NHlg3X5xRXXdsFTvVA+ZSfF1ThJvLtnBMUOzGJSVhMvp4IzROXWeh4FOFrW1ymfr9nDqiF5BJ2qMJo0tKP+i3SHZSU0qqJ22VdQvo+X5i+E51hBpuFDVG1W1IuB4q6qeErYGDS0yqm8aU4dl8fw326gOY2QWQ/swCqoNuJwOqmtr2XPA+jHMjqKCGj8gg9+fOxqwlKb/B7p3mofM5DiKyqqo8tWypbCMwwc0HNZKsOeXKqvrO2RBiZfyqhpG9InesGVbCJyDAktROR1CZnJ8k158O4u9iFjPpyXG9EtjVX4JawvCE6VDRD4RkY8bb2FpzBAUlx+dS0GJl/dXmqnAWMMoqDYQ53RQ7VN22xlto6mgAFxOy9KprqmloMRSUDlpHrKS4ykqq2T7vnJqlYOGIv0OEN6AIT7/EFk0hy3bgt85pThAQaUluEmMc+INsKCKy6tRVfL3V9AzJR63s+V/ef/w5q2vLQuT5PwauNne7gSWAAvC1ZihdY4/tCcDeiQy66st0RbF0AijoNqAyyn4amvrUq5Hc4gPwG0rKF+tkl/sxSGW0kxPdLOvrJrNdgba3EZKJ95lO3wEWFD+dOqDs5IjIXqHcTqEFI+LkkYKyuNyUm5bUMvzihl77we8u7yAwgOVzbqXB3LsIdkMz0nFVxue4R5VXRiwfamqNwHHh6UxQ1A4HcKlkwcyf8s+lucVR1scQwBGQbUBl8NBdY3WKajsIH7wwi0P2BZUcQVZyZaFkBLvoqzKV28VZbZuQW0qLCMxzkmv1Ogq3bZQ6vXx7Fdb8FbXUFxRTWqCG4/bgddXw/7yKp78fCMAK3YWU1RWVRciqjUG9kg8aGFnqBCRHgFbloicBrQ95IchpJw/oT8ZiW7ueGuFmYuKIaKioOzO+aGIrLc/m0zRKiJzRGS/iMyOtIxN4XYKvhrLgopzOUj1RHcZmX+4qrqmlvxiLzn2/Eqyx4WqtZg4xeOq83jzE+8+2ILaXFjGoKykTrly/quNhZR4fZYF5XZSWV3Lj/7xHf9bZoVC6pPmoehAFZlJwSnfOJcjbAoKWIg1pLcQ+Br4FXBluBozBEdagpvffW8US7fv5/Y3lpvI6jFCtCyo24C5qjoMmGsfN8WDwI8iJlUruJ0OfLXWHFR2cnzUf8zrhvhqlF0l3joHgOR4SyGtyi+hf0biQXL6h/gC52r8Cqoz4VfI/nA3fgXlra5h+Y76oZpKXy2FByrJCtKCinc56ta7hRpVHaSqg+3PYap6qqrOC0tjhjZx1pg+/OLEobyyII/fv7vaKKkYIFomwDnUj7vPAj4Fbm1cSFXnisjxjc9HC5dTqK6pZX95NT1jYCjMZVtQvlrLgppsr4ZPti271fklnN5EVG//EJ//R7jKV8v2veWcE+PhjRrz/FWTOOmhz/DVqD0H5SLe7cDbSLkUHqii0lfbYoinQOLdobegROT7LV1X1ddD2qChXfzylEPYX1HN37/YTKrHzfUnNR301BAZoqWgeqmqPxR1AdChlK0icg1wDRDWMPJuhwOfPQc1MDNySQqbl8eyjIorrACZvdOsNT4p8fV/1v49Dl7309iC2rbX9vbL7lwWlF/R+mpr65wkXI565XLbtMN48P21dUFkM5ODHOJzOsNhQZ3dwjUFjIKKAfzZfveXV/Po3PVMG53D0J6dw3GoK9KqghKR3wH3qKrPPk4FHlXVFgNcishHQFP5Gm4PPFBVFZEO2dKq+jTwNMCECRPCZpf7Lag9ByqZkNvktFlE8VtQefusdZ+Bc1B++vc4WJE2tqA6mwefH5etoEsqfNTUKmkJbgLnt3PSPCS4neywF+kG6yThdknILajW+oshdhAR7jp7BB+v2c1DH6zlb5ccEW2Rui3BzEG5gG9FZIyInALMx5rgbRFVPVlVRzWxvQXsEpEcAPtzd0e+RKRwOR1UVNewt6wqKJfl8Mtj/UBvs9NE+Oeg/ItYoWkF1diC8nv7NXZHj3X8ES/8cQfTEtx13w2sqBoet6NOgWcF6SThcgg1IZ5/EJFnA/YvC2nlhpCTmRzPVVMH8d6KApbl7W/9BkNYaFVBqepM4BbgW6z5ojNV9fEOtvs24O+klwFvdbC+iBDnlLqIDdFepAvWwmGot6B628FOewUozxFNJFRsbEFt21tORqK7gWLrDPgtqL1lltt/qsfdIAp7TloCHrezbllAsBaUU4Sa2pAb4oH5nm4IdeUicr6IrBSRWhGZ0OjaTBHZICJrbbd2QxBcecwgMhLdPPj+WuMwESVaVVAicizwGHAvljPDX0Sko7PpDwCniMh64GT7GBGZICLPBLT9BfBf4CQRyYt253I5HOwstpRBLCgovwXlD27qt6BSE+qH+JpaTOxXbH4LKm9fBf0yoj+n1lb8FpQ/pmCyx1UXdR6gZ2p8XVgnIGgnCae9vqw2tEoq3L9wK4DvA58HnhSREcCFwEjgdOCvIhLbuVRihBSPm5+fOIwv1hfy6sK8aIvTLQnGSeJPwPmqugrqvJE+Bg5rb6OqWgSc1MT5BcBVAcdT29tGOHA5Ba+9dijaUSSgfqHudtsC8lsPIsK/Lj+SeLejSVd4h0OIC3ClzttXziG9Ip94saP414H5U4EnxjkbWFAed/1xSrwr6BxX/mhIvlolLnSBc/vZSQolYL8OVf1FRypX1dVAU3/vc4CXVLUS2CwiG4CJWGuwDK1w+dG5fLiqgN++tZLR/dKimmKnOxLMHNRkv3KCOnfYKeETKXYJjOMWCxaUfx1USYAHn58TDuvJ0UOymr033uXAW12DqpK3r6LJuapYx29BlVZaFlSC20VCnKWE/Fai34IKdnjPqte2oEI7rHMz9Yt0/fuBW7joC2wPOM6zzzVARK4RkQUissBkkq3H6RAevfBwUhNcXPGv+ewqCV8yS8PBBGNB3d7MgtR7QyxLzOMKeJvOCtJlOZy4nIEOAW2TJ95luVLvOVBJpa+21TQUsYjT/r88EGBBDbAV7WDbZd5jK6xgh/egoQUVKlR1VkfraMkz1nY+ajeR8oTtjPRK9fDPy4/kgie/5op/zeeVayeTHG+SkUeCYJ5yWcC+BzgLWB0ecWIbt+0hlpHoJs4VrSAc9bgDFGZjC6o1PG4HldU1dQ4WnVFB+XN0BQ7xZafEM3PaYUwfZ02TJthzUm15ofBbUGFwlOgQqnpyO27bAfQPOO5nnzO0gZF90nj84vFcNWsBN/xnMU9fOqFT5E3r7ATjxfdQwHYfVgSIwWGXLAaJt1+tY2F4DxoOOea0kueoMf5wPvUKqvMN8YFlRZZ4rSE+T5wTEeEnxw0hx1bY/nmntnw/e+Q05hRUO3kbuFBE4kVkEDAM+C7KMnVKTji0J3efPYK5a3bz4Ptroy1Ot6A9ZkAi1ltYtyPJNuvbaq2EC78XH9S7mAeLP2adP8pCYAr5zoTLIVTXWIoksQkniLJKy1OxqYgazeF0xqYF1RIicq6I5AGTgf+JyPsAqroSeAVYBcwBrlPVgzM6GoLiR5NzuXjSAJ78bCNvLDaefeEmmEgSy6l3kXUC2XTD+SeoV1A5bVQG4SJwmLG1TLGN8VtQ2/dWkJkUV/fdOhv+YZY4p6PBnJyfkX1S+Wj1LgZnBx8lwz+3FQ4FZVsx1wO5BPQ/VZ3ekXpV9Q3gjWau3Qfc15H6DfXcPX0kG/cc4NbXljNhYI9O6WDUWQjmV+msgH0fsMsf9qi74VcIWSnBT7iHk3hX4KLUtltQlT7LguqM809+/I4rfu+9xlx3wlAmDerB5CGZba4z1NEkbN4E/gG8A5jEQ50Qt9PBIzPGcdyDn/LkZxu579zR0Rapy9KsghKRHvZuaaNLqSKCqu4Nn1ixiT+MTo8gQ+ZEkvZYUKVeH0VlVRzWu/OtgfLjt5oSm1FQcS4HRw9t3t2+KRx+BVUTFgXlVdXHWi9miGVy0hI4c3QOs5flc9fZI2PCaaor0pIFtRBraK8pVxWlGzpKXDElF1+N8sOJ4YuY3lbG9k9nV7GXFE/bwhT5Lai9ZVUx4TLfXlqzoDpSZ5gsqEdF5C7gA6DSf1JVF4WjMUP4OHN0Dm8s3sFXGws5/tCe0RanS9KsglLVQZEUpDOQGOfihpNjKz/MC1dNavINojXiXQ7Kq6xU6emdLAZfIP45qOYsqPZQZ0HVhmUEbjRWEs4TqR/iU/vY0MwtKycAAB2kSURBVIk4ZlgWSXFO5qwoMAoqTDRrl4rIzwP2R0ZGHENbSY53tcvBId7lpPBAJaqQlhgbc2rtwW/tJLpD5+RRZ0GFZ4bofGCwqh6nqifYm1FOnRCP28mJw3vxwapdncrjszPR0sDpjwP2nwu3IIbI4nE76uIKdgULyhNKC8r24vOFx4JaAaSHo2JD5Jk2qjd7y6r4bnO3m5KPCMG+dpol012M+IA1Q+mJnVdB+RcrN7UGqr34lV549BPpwBoRmU/DOagOuZkbosPxh2bjcTuYvWxnmzxFDcHRkoJKF5FzsaysVDuKeR120FhDJ8UT4HXUmRVUOOag/EN8YbKg7gpHpYbokBjn4szRfXhj8Q5uOf2wTpdTLdZpaYjvM2A61jqoz4GzA7azWrjP0AmIb5AnyXjxBeJ3kghxNHMAVPUzYA2QYm+r7XOGTsoVU3Ipr6rhlfnbWy9saBMtefFdEUlBDJElMDV6j07sJBFWCyoM66BE5ALgQazkn4KVAPRmVX015I0ZIsKovmlMHNSDv3+xiYuPGkBiXOeMyhKLmNVl3ZQUj6vJ/c6Gs86CCt138DtJhGkd1O3Akap6mapeipU88M5wNGSIHLeefii7Syt58rNN0RalS2EUVDclLaHeanJ04rQBdYFiQ2lBOcMXiw9wqOrugOMiTD/s9BwxsAdnjcnhyc82smpnSbTF6TKYjtFNyejEjhGBVNlp60O6UDeMwWKBOSLyvohcLiKXA/8D3gtHQ4bIcvf0kWQkurnuxUXsL6+KtjhdglYVlIg4RWS6iPxCRG7yb5EQzhA+Dullxd+bNKhHKyVjm2p7NW1CCN3M6xfqhsVJ4mbgKWCMvT2tqreEvCFDxMlKjucvF41nx74Krpy1gIoqk9WkowRjQb0DXA5kUu951KHooiLSQ0Q+FJH19mdGE2XGicjXIrJSRJaJyIyOtGloSEZSHO/+YipPXzoh2qJ0iDoFFUILyhmgoFbsKGby/XP5emNRSOoWkT+o6uuqepO9vSEifwhJ5YaoM3FQD/584TgWbdvHFc9+R3FFdbRF6tQEo6D6qer3VfUuVb3Hv3Ww3duAuao6DJhrHzemHLhUVUcCpwN/FhGzAj+EjOiT2unXbYRjiM8Z4GZe4q0mv9iLhG6a7pQmzk0LWe2GqHPG6BweuWAcC7fu47y/fcXmwrJoi9RpCUZBvScip4a43XOAWfb+LOB7jQuo6jpVXW/v7wR2YyVLNBjqqKob4gudF5+zbqGuUm5n5O2oAhSRn9rJPw+zRwT822ZgeQdFNsQY3zu8L//+8SR2l1Yy7dHPeeqzjVT6zJBfWwlGQX0DvCEiFSJSIiKlItJRN5Veqppv7xcAvVoqLCITgThgYzPXrxGRBSKyYM+ePR0UzdCZqAyjBVVTq5RX+xVUhxXgi1iL3N+i4aL3I1T14o5Wbog9Jg/J5P0bj+WYoVnc/94aTnroM95cvMMElm0DwSioh4HJQKKqpqpqiqqmtnaTiHwkIiua2M4JLKeqSn1K+abqycEKVnuFqjYZe0ZVn1bVCao6ITvbGFndCf8cVEgVVIAXX3mlLyT1q2qxqm4BfKq6NWDbKyImGHMXpXeah2cuO5LnrpxIqsfNjS8v4cSHPuXZLzdTdKCy9Qq6OcG8Fm4HVtiKJGhU9eTmronILhHJUdV8WwHtbqZcKpYb7u2q+k1b2jd0D/wR2RPbkXKkORpYUFWhGeILoEHqGhFxAUeEqnJDbDJ1WDZTrs/ivRUFPDNvE3e/s4p7Z6/iyNweTBqcyYSBGQzPSSUrOQ4J4YRnZyeYXr0J+FRE3qNh9OWHO9Du28BlwAP251uNC4hIHPAG8G8TBsbQGqEM1xSooCrsIb6OegmKyEzgN0CCPUTu/xWqAp7uUOWGToHDIZw5Joczx+SwYkcx768sYO7q3Tz+8Xr8o34p8S4GZycxMDOJnHQPOakeeqd56J2WQE6ah8ykOFzO7rN8NRgFtdne4uwtFDwAvCIiVwJbgQsARGQCcK2qXmWfOxbItBc0AlyuqktCJIOhC+DPaxWulO/lVT4cAnEd/FFQ1fuB+0XkflWdGQIxDZ2YUX3TGNU3jV+deiil3mqWbi9m/e5SNheWsbmwjMXb9zFnRWWdE1AgiXFOkuNdpHhcJHvcpHpcJMe7SPW4SU90k5bopl9GIsN7pzAoK6lTK7RWFZTfpVxEku3jAx1tVFWLgJOaOL8AuMrefx54vqNtGbo2791wLBt2d/hfsgGOQAuqqpbEOFfIhl1UdaaITMd6+QL4VFVnh6RyQ6ckxePmmGFZHDMsq8F5VWVvWRX5xV4Kir0UlHgpPFDJAa+PA5U+Sr0+Sit9HPBWU1DspcRbzb7y6rqlF2C9WA3pmczw3ikclpPC6SNzGJCZGOmv2G5aVVAiMgrLSaGHfVyItT5pZZhlMxhaZVBWEoOykkJap6vBEJ8vpNaZiNyPFSD2BfvUDSJytKr+JmSNGLoEIkJmcjyZyfGM6psW9H0VVTVsLixj7a4S1uSXsqaglK82FvH64h088N4arjxmEDedcmhI/6/DRTBDfE8DN6nqJwAicjzwd+DoMMplMESNhhZUTUjDKAFnAuP8HqkiMgtYjDU/ZTB0mIQ4JyP6pDKiTyocXn++oNjLo3PX8fcvNvP+yl3c//3RTBma1XxFMUAwg5NJfuUEoKqfAqF9ZTUYYghXIy++ULqw2wRGRAn+1dhg6AC90zzc//0xvHTNUTgELn7mW+54czltdNCOKMEoqE0icqeI5NrbHViefQZDl8Qfzdxne/F5QmtB3Q8sFpFnbetpIXBfRysVkQdFZI0dneKNwLBgIjJTRDaIyFoROa2jbRk6N0cNzmTOjcdyxZRcnv9mG7O+2hJtkZolGAX1Y6wQQ68DrwFZgMm2a+iy1MXis4f4QmlBqep/gKOo70+TVfXlEFT9ITBKVccA64CZACIyArgQa/3V6cBfRST2Jx8MYcXjdvLbs0YwdVgWD3+4jn1lsZkeJBgFdbKq/kJVx6vqEap6I00HvDQYugTORhZUqOagRMQlImKH+VqMtWwjJxR1q+oHquqzD78B+tn75wAvqWqlqm4GNmA5aRi6OSLCnWeN4EClj39+uTna4jRJMAqqqTUbZh2HocvicAgiVjTziqqakHg7icjVWBFTttr7c4HzgJdE5NYON9CQH1OfBLEvVjQYP3n2ucbymXiW3ZBDeqVw9JAs3l2e33rhKNCsF5+ITAPOAPqKyGMBl1IBX9N3GQxdA5dDrGjmofPiuxEYgpVLbTUwUFULRSQRmA+0mhNKRD4Cejdx6XZVfcsucztW/3yhiXLNoqpPY0e0mDBhQuzOmhtCzmkje3HnWyvZsLuUoT07lOov5LRkQe0EFgBerIlc//Y2YCZaDV0ah4g1B1UdsjmoKlXdp6rbgA2qWgigquVY4Y5aRVVPVtVRTWx+5XQ5cBZwcUDszB1A/4Bq+tnnDAYAThlhvfPMWVEQZUkOplkLSlWXAktFpJeqzgq8JiI3AI+GWziDIVrEuxx4q2vsIb6QBKJNEJHDsV4K4+x9sTdPRysXkdOBW4DjbKXn523gRRF5GOgDDAO+62h7hq5D7zQPhw9I5/2Vu/j5icOiLU4DgpmDurCJc5eHWA6DIabISIqj8EAVVTW1oRriy8dKXfMnrBxoDwMPBRx3lMexhg8/FJElIvIkgB3x5RVgFTAHuE5VTeY8QwOmjerN8h3FLN2+P9qiNKClOaiLgB8Cg0Tk7YBLqcDecAtmMESTjMQ4duyvAEKTakNVT+hwJS3XP7SFa/cRgrVWhq7LRRMH8PTnm7nzrRW88pPJoV77125aGrv4CuutLwvrTc9PKbAsnEIZDNGmR1IcK3YUA+DpBDHLDIaOkOJx83/fG8m1zy/i5leX8ecZ4+rWA0aTZof47Gyfn6rqZGAN1vBBCpAXsN7CYOiSZCTGsbvUSn+WGCNvkwZDODl9VA63TTuMd5bu5HezV0VbHCCIOSgROR9rUvV8rBxN34rIeeEWzGCIJime+sGFzhD12WAIBdceN4RLjhrAv7/ewrai8lbLh5tgnCTuAI5U1ctU9VKsVeh3hlcsgyG6BA5veNydN+GbwdBWrj9xGC6Hgyc/3xhtUYJSUA5V3R1wXBTkfQZDpyVQQbnDnJFURBaFtQGDoQ30SvVw3oR+/HfBdnbajkLRIpieN0dE3heRy+2FgP8D3g2vWAZDdHEEZNB1OcKroFR1fFgbMBjayM+OHwLAY3PXR1WOVnueqt4MPAWMsbenVTXUscMMhpgi0IEpzhV6byYRSRWRHv4t5A0YDB2gX0Yilxw1kFcWbGdNQUnU5Ajq1VBVX1fVm4DfA2+GVySDIfoEDvGF0oISkZ+ISAHWUg1/+LAFIWvAYAgRN5w0zHI/n706akkNm+15InKUiHwqIq+LyOEisgJYAeyyw6q0G/ut8UMRWW9/ZjRRZqCILLJXxa8UkWs70qbB0BYaDPE5Q2pB/Rorb1Ouqg6yt8GhbMBgCAXpiXHccNIw5m0o5NO10Ylw39Kr4eNYFtN/gI+Bq1S1N3AsVlbQjnAbMFdVh2GlHbitiTL5WMncxgGTgNtEpE8H2zUYgiJQQcWF1kliIxB9/12DIQguOWogg7KSuP+96FhRLfU8l50E7b9Agap+A6Cqa0LQ7jmAPwDtLOB7jQuoapWqVtqH8a3IajCElECd5AqtgpoJfCUiT4nIY/4tlA0YDKEizuXg+hOHsm7XAb7cUBTx9lvqebUB+419DTuqSnvZWUXBCpTZq6lCItJfRJZhJVz7g6rubKacSbZmCCmOBm7mIR3iewprROIbGqaxMRhikjPH5JCZFMezX22JeNstxeIbKyIlWOkAEux9CDI9QEvJ1QIPVFVFpEmFp6rbgTH20N6bIvKqqu5qopxJtmYIKYFDfCFeB+W2HY4Mhk5BvMvJDycN4PFPNrA6v4ThOakRa7ulWHxOVU1V1RRVddn7/mN3axW3klxtl4jkANifu1upayeWg8bUtn09g6F9OMOnoN6zLf4c42Zu6Cxcdcxg0hPc3P7Gcqpralu/IUREa17nbeAye/8y4K3GBUSkn4gk2PsZwDHA2ohJaOjWBA7xhdiL7yLseSiMm7mhk5CW6Oaec0axaNt+7nxzBbW1kRmoCkmq0HbwAPCKiFwJbMUKQouITACuVdWrgOHAQ/bwnwB/UtXlUZLX0M1osFA3hBaUqg4KWWUGQwSZPrYP63eV8pePN3Cg0sf93x9NiqfVwbQOERUFpapFwElNnF8AXGXvf4gVucJgiDgNF+qGzoISkUubOq+q/w5ZIwZDmLjplENIinfxxzlrWLR1H3dNH8mpI3ohEp7cUcZ122BogkAniRAnbjsyYJsK3A1MD2UDBkO4EBGuPW4Ir/70aFI8bn7y3EJ+/Ox8thaVhaU9o6AMhiYIVFChfDtU1esDtquB8UByyBowGCLA+AEZzP7FMdxx5nC+27yXMx+bx7z1hSFvxygog6EJwpxhI5AywMxLGTodbqeDq6YO5oObjqNvegKX/+s7Xl+UF9I2ouUkYTDENI4wjamLyDvUL3R3ACOAV8LSmMEQAfqmJ/DKtZO59rmF3PTKUgpKvPzs+KEhqdsoKIOhCcKloIA/Bez7gK2qGtrXToMhwqQluHn2x0dy83+X8cc5a+mRGMeFEwd0uF6joAyGJgixYwQiMhQrxNdnjc5PEZF4VY1+fm2DoQPEu5w8fMFY9pVXcdfbKzlyUA+GZHdserVbKKjq6mry8vLwer3RFiWseDwe+vXrh9sd3rUJ3QFHiBUU8GesBbqNKbGvnR3qBg2GSONyOnjo/LGc8sjn3PrqMl796dEdqy9EcsU0eXl5pKSkkJubGzZ//WijqhQVFZGXl8egQWbOvaOEXj/Rq6mF5qq6XERyQ96awRAleqZ6uOGkYdw7exXrd5UyrFdKu+vqFl58Xq+XzMzMLqucwHKFzszM7PJWYqRwhv5/Jb2FawmhbsxgiCanj7LihH+5oWOu591CQUFo17LEKt3hO0YK/xBfCB/pAhG5uvFJEbmKEKTbEJHficgyOwP1B/7knmLxmIhssK+P72hbBkNr9ElPoGdKPEvzijtUT7cY4jMY2orfiy+E3nw3Am+IyMXUK6QJQBxwbgjqf1BV7wQQkV8AvwWuBaYBw+xtEvA3+9NgCCtj+6ezdPv+DtXRbSyoWGTLli0kJCQwbtw4xo4dy9FHH83atS0HbF+yZAnvvvtuhCTsvvgX6oZqqE9Vd6nq0cA9wBZ7u0dVJ6tqQQjqLwk4TKJ+rdU5wL/V4hsg3Z/qxmAIJ+P6p7OpsIziiup212EUVJQZMmQIS5YsYenSpVx22WX8/ve/b7G8UVCRoc6CCnEPUdVPVPUv9vZxKOsWkftEZDtwMZYFBdAXKyO1nzz7XON7TVZqQ0g5rLflHLFh94F219HthvjueWclq3aWtF6wDYzok8pdZ49sscxtt91G//79ue666wC4++67SU5uuEagpKSEjIwMwHLs+OlPf8qCBQtwuVw8/PDDTJkyhd/+9rdUVFQwb948Zs6cyYwZM0L6XQwWfgUVBmeJdtNSlmpVfUtVbwduF5GZwM+Bu4Kt22SlNoQa/xqojXsOcMTAjHbV0e0UVLSYMWMGN954Y52CeuWVV3jqqae48847GTduHKWlpZSXl/Ptt98C8MQTTyAiLF++nDVr1nDqqaeybt067r33XhYsWMDjjz8eza/T5fEv1A1jRIk2o6onB1n0BeBdLAW1A+gfcK2ffc5gCCv9MhKIczrYuMdYUEHTmqUTLg4//HB2797Nzp072bNnDxkZGfTv379uiA/g5Zdf5pprrmHOnDnMmzeP66+/HoDDDjuMgQMHsm7duqjI3h2pH+KLHQXVEiIyTFXX24fnAGvs/beBn4vIS1jOEcWqmh8NGQ3dC5fTwcDMRDbtaX8qjm6noKLJ+eefz6uvvkpBQUGTQ3PTp0/niiuuiIJkhsb49VKoQx6FkQdE5FCgFitL9bX2+XeBM4ANQDlg/sEMEWNIdjLrdpe2+36joCLIjBkzuPrqqyksLOSzzz6jsrKywfV58+YxZMgQAKZOncoLL7zAiSeeyLp169i2bRuHHnoo69evp7S0/X9wQ3DUD/FFWZAgUdUfNHNegesiLI7BAMCQnkl8uHoXVb5a4lxt9zgyXnwRZOTIkZSWltK3b19ycixP340bN9a5mf/mN7/hmWeeAeBnP/sZtbW1jB49mhkzZvDss88SHx/PCSecwKpVqxg3bhwvv/xyNL9Ol8YRg3NQBkNn45BeKdTUKut2te+l2lhQEWb58vpwbLm5uVRUVDRZzuPx8K9//eug8z169GD+/Plhk89gUefF11lMKIMhBhk/wPLeW7xtH6P6prX5fmNBGQxN4Ax9JAmDodvRLyOBrOR4Fm9rX0SJqCgoEekhIh+KyHr7s1kneRFJFZE8ETF+1YaI4V+gG+qFugZDd0JEGD8gnflb97br/mh1v9uAuao6DJhrHzfH74DPO9qgNVfctekO3zFSuGzN1DfdBBo3GDrClKFZbN9bwZbCtrubR0tBnQPMsvdnAd9rqpCIHAH0Aj7oSGMej4eioqIu/QPuzwfl8XiiLUqXYGjPZC6a2J8bTz4k2qIYDJ2a4w7JZmy/NPaVV7X53mg5SfQKWCxYgKWEGiAiDuAh4BKgxRX0InINcA3AgAEDDrrer18/8vLy6OoxxvwZdQ0dx+kQ7v/+mGiLYTB0enKzknjr58e0696wKaiW4oYFHqiqikhTps3PgHdVNa+1PEetxRFzu90my6zBYDB0MsKmoFqKGyYiu0QkR1Xz7dD/u5soNhmYKiI/A5KBOBE5oKotzVcZDAaDoYsQrSG+t4HLgAfsz7caF1DVi/37InI5MMEoJ4PBYOg+RMtJ4gHgFBFZjzW/9ACAiEwQkWeiJJPBYDAYYgjpap5tIrIHK1hmU2QBhREUpy0Y2dpHLMg2UFWzoyxDyDB9KCwY2VqmyT7U5RRUS4jIAlWdEG05msLI1j5iWbauSCw/byNb+4hl2cw6eYPBYDDEJEZBGQwGgyEm6W4K6uloC9ACRrb2EcuydUVi+Xkb2dpHzMrWreagDAaDwdB56G4WlMFgMBg6CUZBGQwGgyEm6TYKSkROF5G1IrJBRCIekUJE+ovIJyKySkRWisgN9vkmc2OJxWO2vMtEZHyY5XOKyGIRmW0fDxKRb+32XxaROPt8vH28wb6eG0657DbTReRVEVkjIqtFZHKsPLfuhOlDrcoXk32oM/efbqGgRMQJPAFMA0YAF4nIiAiL4QN+paojgKOA62wZmsuNNQ0YZm/XAH8Ls3w3AKsDjv8APKKqQ4F9wJX2+SuBffb5R+xy4eZRYI6qHgaMteWMlefWLTB9KChitQ913v6jql1+wwo8+37A8UxgZpRlegs4BVgL5NjncoC19v5TwEUB5evKhUGWflj/pCcCswHBWlnuavz8gPeByfa+yy4nYXxOacDmxm3EwnPrTpvpQ63KEpN9qLP3n25hQQF9ge0Bx3n2uahgm/SHA9/SfG6sSMr8Z+AWoNY+zgT2q6qvibbr5LKvF9vlw8UgYA/wL3v45BkRSSI2nlt3Iqaeq+lDQdOp+093UVAxg4gkA68BN6pqSeA1tV5ZIur3LyJnAbtVdWEk220DLmA88DdVPRwoo344AojOczNED9OH2kSn7j/dRUHtAPoHHPezz0UUEXFjdawXVPV1+/QusXJiIQ1zY0VK5inAdBHZAryENUTxKJAuIv50LIFt18llX08DisIgl588IE9Vv7WPX8XqcNF+bt2NmHiupg+1mU7df7qLgpoPDLO9auKAC7FyUkUMERHgH8BqVX044JI/NxY0zI31NnCp7VVzFFAcYJKHDFWdqar9VDUX67l8rFYurk+A85qRyy/veXb5sL19qWoBsF1EDrVPnQSsIsrPrRti+lAzxHIf6vT9J1qTX5HegDOAdcBG4PYotH8Mlhm9DFhib2dgjT3PBdYDHwE97PKC5TW1EViOlbAx3DIeD8y29wcD3wEbgP8C8fZ5j328wb4+OAJyjQMW2M/uTSAjlp5bd9lMHwpKxpjrQ525/5hQRwaDwWCISbrLEJ/BYDAYOhlGQRkMBoMhJjEKymAwGAwxiVFQBoPBYIhJjIIyGAwGQ0xiFFQMIiKZIrLE3gpEZEfA8VdhaO9yEdkjIs80c/1TEZkQwvYetL/Xr0NVp8EQiOlDXQNX60UMkUZVi7DWLiAidwMHVPVPYW72ZVX9eZjbAEBVbxaRski0ZeiemD7UNTAWVCdDRA7Yn8eLyGci8paIbBKRB0TkYhH5TkSWi8gQu1y2iLwmIvPtbUoQbSSIyEti5Y55A0gIuPY3EVkgVj6ee+xzJ4rImwFlThGRN8TKj/OsiKywZfplyB+IwdBGTB/qPBgLqnMzFhgO7AU2Ac+o6kSxErldD9yIFRPsEVWdJyIDsEL9D2+l3p8C5ao6XETGAIsCrt2uqnvFyg80177+CfBXEclW1T3AFcA/sd5g+6rqKLASp4XoexsMocL0oRjGWFCdm/mqmq+qlVihST6wzy8Hcu39k4HHRWQJVpytVLGiQbfEscDzAKq6DCtEip8LRGQRsBgYCYxQKxzJc8AldgeaDLyH1eEHi8hfROR0oEHkaYMhBjB9KIYxFlTnpjJgvzbguJb6v60DOEpVvR1tTEQGAb8GjlTVfSLyLFZcMYB/Ae8AXuC/auW52SciY4HTgGuBC4Afd1QOgyGEmD4UwxgLquvzAdZQBQAiMi6Iez4HfmiXHwWMsc+nYuWTKRaRXljpoQFQ1Z3ATuAOrI6GiGQBDlV9zT4/vqNfxmCIAqYPRQljQXV9fgE8ISLLsP7en2O9ibXE37AycK4GVgMLAVR1qYgsBtZgZd38stF9LwDZqrraPu5r1+N/EZrZ0S9jMEQB04eihIlmbkBELscKq98hF1kReRxYrKr/CKLs3UTG9ddgCDumD4UHM8RnAKgApkkziwyDQUQWYg1jPB9E2QeBS7CGOgyGroDpQ2HAWFAGg8FgiEmMBWUwGAyGmMQoKIPBYDDEJEZBGQwGgyEmMQrKYDAYDDGJUVAGg8FgiEn+Hyo7HR+lhCpkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x180 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "axes = ml.plots.water_flow(\"Bottom Flux\", figsize=(6, 2.5))\n",
    "#plt.savefig(\"../../figures/water_flow.eps\", bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9c. Plot the water content over time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[Text(-100.0, 0, 'Jan-07'), Text(0.0, 0, 'Jan-08'), Text(100.0, 0, 'Jan-09')]"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAADCCAYAAAA7BNegAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2df5QlVXXvP7u7p4draGYCQxv50d34A7UnSyRMTGiCPzKYhRORcYUk6hMfYhpNFmElJhMlklk6mERvgkSCeYooSnRpYow88UeQ7kQURp80KCIDMzDg8GMEYXCGGR1+zLDfH3WqOV1TdW/VvVW36lbvz1pndd+quqdO1T11vmfvs88pUVUMwzAMo24MlF0AwzAMwygCEzjDMAyjlpjAGYZhGLXEBM4wDMOoJSZwhmEYRi0xgTMMwzBqiQmcYRiGUUtM4IzSEJHzReTrkW13Jmx7Q4r83isin8m7nHkjIhMioiIylFN+Z4nI9V18/5Uicn/G70yKyJyI/MylGRGZ9PaLiHxQRHa49EERkU7LaBidYAJnlMm3gCkRGQQQkecAS4DjI9ue744tlLwEZ5GwHTgDOBRYAXwZ+Ly3/xxgLXAc8BLgNODtPS6jscgxgTPK5EYCQXup+3wy8D/A5si2raq6HUBEPiwi94nIYyJyk4ic7LafCvw18IciskdEbnHbl4nIJ0TkJyLygIi83xPPs0TkBhG5WER2AO+NFlBEBkXkr0Vkq4jsduc82u2bEpEbRWSX+zvlfe+bInKhy3+3iHxDRFa43aFY73RlPdF952wRud1ZRNeIyLiXn4rIO5w1u1NEPuKspBcDHwVOdHntjLvRIvJWl/duEblbRN7utv8S8HXgCPf9PSJyRLsfTlV3quqPNVgKSYD9BB2RkP8NXKSq96vqA8BFwFlJ+YnI6SLyA/e7bnW/Z3gf3y8iG13ZrhaRw0Tks+7YG0Vkol15jUWKqlqyVFoiELQ/d/9fCpwN/G1k2ye9498MHAYMAX8BPAgc5Pa9F/hMJP8vAR8DfgkYBb4HvN3tOwvYB/ypy68RU751wK3ACwka8uPc+Q8Ffgac6b77Rvf5MPe9bwJbgWOBhvv8AbdvAlBgyDvP6cBdwItdfhcAG739CnwFWA6MAQ8Dp3rXcX2b+/y7wPPcNbwC+AXwa27fK4H7O/z9drp7+DRwgbd9F/Ab3udVwO6EPF7mjn81Qaf7SOBF3n28y5V9GbAJ2AKc4u7TlcAVZddjS9VMZsEZZXMd8HL3/8nAt13yt10XHqyqn1HVHaq6T1UvApYSiM8BiMizgTXAn6nqz1X1p8DFgD+et11V/9nltzcmmz8iaLg3a8AtqrqDQDDuVNV/dd/9HHAHgSsu5ApV3eLy/XeesUrjeAfw96p6u6ruA/4OeKlvxREI5E5VvZegY9AqvwWo6ldVdau7huuAbxDc265Q1eUEwnMu8H1v18EEohWyCzg4YRzubQSdmGtV9WlVfUBV7/D2X+HKvovA2tyqqjPuPn0BOL7b6zDqiQmcUTbfAn5LRA4FDlfVO4GNBGNzhwK/ijf+JiJ/6Vxtu5w7bhnBGFAc4wQu0J84t95OAmtu1DvmvjblO5rAEotyBLAtsm0bgfUR8qD3/y8IGv0kxoEPe+V8lMDa6jS/BYjIa0TkuyLyqMt/Dcn3LfrdMc99uSe6X1V/TuAmvVJEwnu7BzjEO+wQYI+qxq3unnSPQx7y/t8b8zn1fTAWFyZwRtl8h0CkpoEbAFT1MYIghmkCC+seADfe9lfAHwC/7KyHXQRCAIEbz+c+4Alghaoud+kQVV3pHdPudRr3EbjHomwnECWfMeCBNvklnfM+Atfpci81VHVjh/nNIyJLgS8C/wg82923r5F83xZmrnqvqh4cpoTDBoBn8Ywg30bgzg05zm2LI+keG0ZXmMAZpeLcd3PAOwlckyHXu21+9OQIwXjPw8CQiKxnoZXwEDAhIgMu758QuOIuEpFDRGRARJ4nIq/IUMTLgQtF5AUuqOMlInIYgUAcKyJvEpEhEflDYJJgnKwdDxOMWT3X2/ZR4HwRWQnzwTG/n7KMDwFHichwwv5hAlfuw8A+EXkN8DuR7x8mIstSng8RebWIHO+CcA4BPkQwBnm7O+RK4J0icqQLWvkL4FMJ2X0CeKuIrHa/0ZEi8qK0ZTGMJEzgjCpwHYHb0J/L9W23zRe4a4D/Iggy2AY8zkIX4xfc3x0icrP7/y0EDfwmggb4P4DnZCjbhwjGz74BPEbQGDfcONxrCRruHQSW5WtV9ZF2GarqLwgCaW5wLsnfVNUvAR8EPi8ijwE/Al6Tsoz/TWAdPSgiB5xfVXcD57nr+BnwJoKw/nD/HcDngLtdedpGURIEu3yOwILeSmCBnaqqj7v9HwOuJgjQ+RHwVbftAFT1e8BbCcZHdxHUh6h1bBiZkXiXuGEYhmH0N2bBGYZhGLXEBM4wDMOoJSZwhmEYRi0xgTMMwzBqSd8vLrtixQqdmJgouxiGYRhGh9x0002PqOrheefb9wI3MTHB3Nxc2cUwDMMwOkREoqsC5YK5KA3DMEpidnaWk046idnZ2bKLUjhlXGvfz4NbtWqVmgVnGEY/ctJJJ7Fx40ampqa44YYbyi5OoaxcuZJNmzYxOTnJbbctXLVNRG5S1VV5n9MsOKPWLKYestF/bNiwgampKTZs2FB2UQqh7OfPLDij1iymHrJhVInZ2Vle//rXs3v37nkRP++88wC45JJLWL169fyxZsH1MWX3YhYzGzZsYHJykp07d9r9NypHnduG9evXs3v3bgYGBli7di2rV69m+fLlbNq0ifXr1/ekDCZwPeC8885j48aN870Xo3vSNgxlPFSGkYbQwtm4cWMt6+aGDRsYGRnh6aef5qqrrgJg7dq1jIyMsHbt2t4UouxXinebTjjhBK06k5OTCujk5GTZRakNU1NTCujU1FTbY2dmZnRqakpnZmZ6UDLDSEdYh0dGRnKtm1Wq7zMzMzo5OaljY2M6OTk53xZGrxmY0wL0oXSB6jb1g8BVqcLVhWazqSMjI9psNssuSkeED/7k5KTVi0VKlnYhy7FZOn+9ICwPoGNjYzoyMnJA+WohcMCpwGbgLuDdMfvfSfDerh8Cs8B4uzz7QeBUTeTypp8tuGazqQMDA/MPfVUaIqN44jpmaepnaPk0Go229bhqnb+ZmRltNBrzXqy48vW9wAGDBC9GfC7BCyhvASYjx7wKeJb7/4+Bf2uXb5UFzq+4VetV9TtJjULc9qrd+7AHGzZYVWmIjPyJ1sfwtx8ZGZk/Jk39DAUuzVBH1eq7anxbODAwMF/36yBwJwLXeJ/PB85vcfzxwA3t8q2ywPljb1XrVdWVuIe7au7AZrOpjUZDly5dqoCKiI6Pj1eibEa+ROtjpxZc1ApqRdXqe5SZmZl5D8bAwIDOzMzUQuDOAC73Pp8JXNri+EuBCxL2nQPMAXNjY2P53PUC8AWuir2qfieuYUhqLDq9/0W5N/1xibQ9c6P/SCteaepYFuGqenvTbDZVROa9GMBmXSwCB7wZ+C6wtF2+Vbbg/Ipb9V5VP5JkrcU1FlksaD/ya3h4OFfxCcvXbDZ1ampKp6en5x90ETELv4ZE62T0cxYxandsmPf09HTlPUa+2xXYq30ucKlclMApwO3AaJp8qyxwqjYOVyRxDUdchJZqZ42In9IM7qchSZSjLhujPkR/8+jnLJ3faEctSSyTnoOi6KQD77tdgf3a5wI3BNwNHOMFmayMHHO8C0R5Qdp8qy5wNg5XHEkPd6PROOBBy/IARh48XbJkSW5WXCsL03fZmMhVn5mZGR0bG2sbKNTOglNt3QFr1UmOzqXzPQS99Bj5ncLwOUnrng3FWPtZ4DQQsDXAFidi73HbNgCvc//PAA8BP3Dpy+3y7CeBK2pi52Ilric8NTU1f8+jjUXaCfehGIYi5/8t8nfzXTYmctXH/738qMgoaQSuVefXr+fRjlqr4JNeLjDhd9BCV3tar4kLMtmj/S5wRaSqC1x0HC5tJJTRniSrLKmx8BukVlGL4YMZdkrCqMcsv1sOLhsbk6s4WetTkosyaVtIVBCjwpU0Ry7L3Lluibr1wykAaQO06PcoyqJS1QUuSi8r3WIgS2Phj3W1Equ4Hnaa3y3OlZS1MxMVOcBErqK08xhEj/MtuGjnJ0uHyK+LoSsyrDO+d6iXHWo/uKUTV7sJXE0Ezqy47vEbjKR5RUmNRadjXf7vFve96P5odOT09HSmZZnGxsYWiJzNk6suacbW21lhSduS8grrWljHfJGLy7OTDnV4XdPT0y3HGpOuLUsbV5rAAYemSMuLKFya1G8Cp/pMBTAXVHaikZJJ1lqrxqLTTkYkrFkHBwd16dKlC8YbfPei/7kTsWo2mx1/13iGIuYyZo2Ojh7TSuDSiFHceG1cnt10qP0Vd3zXY7RscWPhWc9ZpsA9ThD9eE+LdG8RhUuTsgpcEZGMWcdbLCy8c+KixuLufbvecCedjDj3YZjWrFmjjUZDh4aGFgQeTE9Pxx6fNnzbtzhtbK4zssyXTItfv9I8/3Hh/XFuyiyrlUSPbfcsZPVYRL0IUUFNurbw+1nub5kC9/08jikqZRW4sFfSTliy/ECdjLf4C+5aVGV6Wo2PRXuu/ms64lyKcZ2Mdr97+OAPDQ3p4ODgAQ9/2MD45wyDVJYvX67Dw8OZrbDwnFGhGx0dtYUDUpA1LD8N0TrXLr+4COq472SJfEx7XZ1YVP6Um/Dv+Pj4gk7WmjVrdGRkZH57N3PuyhS4g2K2/Uq7Y3qVOrHgfL91Ell6+O3GZ9J8z8bjFpKlg5HUKPhunLiHL64uZG34fDdi0dZ4q171kiVLzH3ZAj8IIrSsO50nFmcptfMMJVlc0TreSfRt3HVGv5s24jOajx8J2cqD0W3dK03gYr8ENxdRmE5SJ2NwaUx2v0KkabiiFThtA23jcfFkEZpuXDP+/W83mJ5Eryfwx7ktzX3ZmugY6cjISMdWXNbIXX/+WprObBHz13yPRdaOeDSfuE5Wq3mAaaiawJXmkoymLALn90raWU7RHzLLwG+4kkaah8fG4+KJila7DkOSmzLNSgr+g591jKwswvszOjqqQ0NDBwhetEcd1ufF+nqeaKdgdHRUR0dHO7ofWSyv6CTtLJ3evKcSRe9BN/mHnbrQTdltnaqawP1JEYXpJGUROL+ypbWcsoS8+lbc2NhYahdbWrfpYiIaLRn+Dkn3tdtwaF8c+nGOYlzPOqzbUdeSiCy6sTvfgou62bJ0ZtKuOBLSSTBLkUMXcXXBv5aw89/r8d1KCVyVUicWXGgR+D33NL7zsDI0m00dHh6OHffotKG1CeALCRuLcEUEv0OS5ArqZoHZ8MHu93GsOPfl6OhorDtzsdS1aMDR9PT0fCOeFISURNwLS0PydFUX2R7EjaWNj48f8Kb5XnoyShO4NONtZY7JdTMPLilAINrr8sVQRA5oQHyrotPeVxkBJ530LnuFL1ih6yd8dU3SQ1/l6+klSeMkoTtusS0H1mqprE6CipI6QnmuNdtp4FqW/JMClsroAJUpcHuBH7ZIt9JH8+CixK1sEVfpk+YjxfV2+sWKK2IwO0/Cnnd00eOqlrdqtOrAxfXg69oxiOuwtloJpx2tAkryrKNxk7nzJq5dyzK8khdlCtx4inRUEYVLk/JYySS6QkWS/zk6HpcUINKpu6yXVly/TFHwe8VZFm81ApIa8Lge/GJxWfp0YnW18hT4ndRu62u0I1KUuzCsC53M0cwLG4MrUODierRJ6w2mnf/SqbusV1acP5GzygEH9g69Yskzsq4fydsVGNeWdNOBTIrOrBsmcAnp2GOPnV9ZopvJrmEvxn/YywgV75VlFQpw2IuvqhWX57iGEU+0UV5sIpclUjprflV+tqqECVxCSprw2o3Qld1j6uXk734Yh+smOtJIR7vw8bT0YzRqEVZc2S6/fqPMMbgJ4B+A/wQuB84FxosoTEcX0CLoo18jxJLG9oo6V9mC3g6LjuwNce61LM9PdOpNP1mCRUctGq0pU+BuAd4GrAa2AR93fz8CLC2iUFnSQQcdtGC1i1aL0vZTkEIvJ393a8XZOFl9iAs+SWuFRJfD6pVY5NUByttVaaSnTIH7kff/993fIWAd8OkiCpUlxQWZtJvj0S+url4FnHR7nrRvaDD6h7h30Q0ODi54B17cd8KFjH1LsOodtBCz4sqjTIH7J+Bc9//NkX1biihUltQqijIcD/Df0dVPjXARASdJ727q9DxRq7lfOg9Ge1ot6gwHvsEguu5idAWgIigyCtKsuN5RpsANAO8B5oAHgXOANzsX5XcynQxOBTYDdwHvjtn/cuBmYB9wRpo800wTCBvhflxothPrqtXUhaTlyTq14sJGbWxszNyUNSR8dsL338W9Ay90Y0aHAKL1rYiAi7D+5enOt2Xzek9pAjd/IDwLeB2wHrgY+GNgeYbvDwJbgecCw25sbzJyzATwEuDKPAWun+mkR+m/at4XnOgYSdLqFlmXGPMXQ7Zw/voTrsWaJqgr6urM25oL61/40s28rTgTud5QusB1fSI4EbjG+3w+cH7CsZ8ygXuGLD3KuPHHsEEJG4Pp6en5nnXc24Q7eagtnH/xERfUlbQIcVQM87b0i1wmy1yVxVMHgTsDuNz7fCZwacKxLQXOuUnngLmxsbF87nCFySIevsuwVYMSF3nWrUhZNOXiJM0QQFzHa3BwMNc3kefpWuzVGKIRYAK3iC041fTiEV1INmuvuRuRSrPqSD+PhxrdExedGaak6My05G3FtRqzNvKlDgJnLsouSPtaj6hlFm1Q2vWWs74+JHruVg1MdCJxUkSrTeyuN+EY3uDgYOw4XjeLDuS9Mk/SGxmMfKmDwA0BdwPH8EyQycqEY03gIrRacSSuwYi+6sdvRFqJV7fi0spNlHYicDcia/QXfpRmHos+F7Eyj/8iUKuTxdD3AhdcA2uALQTRlO9x2zYAr3P//zpwP/BzYAdwW7s8F4vAqcY3/HEun7ieZpY1ArsRuVZWnF+GMOotbozDLLjFSV6LPkfn4+VRl8y1Xiy1ELgi0mISuOjDGh0jaLXKRBa6taCS3ERJE4HTuk+N+hOtF524Gv3nJE9vQBiAFRcpanSHCZwJ3AKKjPLqtteb5CaKE+i4JdWSXjhrLA7yDBbJ0xtgUcLFYQJnArfgYS1ynk4ejUKWkO2kJaFsvGPx0s1Ymrm4+w8TOBO4eXfLwMCAjo6O5jbnJ+k83QhMVgGOThq2FSSMTuuhBSn1HyZwJnA9WdsvPM/k5KSOjo62HVRP6i13utyR9b6NED86OGlsOWnBgnCb1af+wATOBE5VF0ZNFjnY7Yf0DwwM6PT0tA4PD883NNH1CON6y7Z8l9EN0WklcfW9lbVm4f39gwmcCdw8vRjsTgoAiUutJsD661/aAL2RBX+OXNKSXkkWmk3Q7i9M4EzgSsEfS4tL7aYlRIXSQqyNTvDrYSh2rTpNSW/UMKqJCZwJXClEQ7bDhibt1ISom8kaG6MTWnW0kt5gYB6D/sEEzgSuNLoZtLcVIIw8SFrSyzpN9aAogZMg7/5l1apVOjc3V3YxDMPoEbOzs5x99tk8/PDDvO9972PdunVlF8noEhG5SVVX5Z3vUN4ZGoZhFMnq1avZtm1b2cUw+oC+t+BEZDewuexy5MQyYFfZhciBulwHVONaVgCPdJlHFa4jL+xaqke31/FCVR3JqzAhdbDgNhdh2paBiFymqueUXY5uqct1QDWuRUTmuq3jVbiOvLBrqR7dXoeIFDLONFBEpkbHXF12AXKiLtcB9bmWulwH2LVUkUpeRx1clF33bg2jylgdN+pOUXW8DhbcZWUXwDAKxuq4UXcKqeN9b8EZhmEYRhx1sOAMwzAM4wBM4AzDMIxaYgJnGIZh1BITOMMwDKOWmMAZhmEYtcQEzjAMw6glJnCGYRhGLTGBMwzDMGqJCZxhGIZRS0zgDMMwjFpiAmeUhoicLyJfj2y7M2HbG1Lk914R+Uze5cwbEZkQERWRXF5XJSJnicj1XXz/lSJyfwffe5aI/IuIPCIiu0TkW94+EZEPisgOlz4oItJpGQ2jE+rwPjijf/kW8G4RGVTV/SLyHGAJcHxk2/PdsYUiIkOquq/o89SIywjakBcDjwIv9fadA6wFjgMUuBa4B/hoj8toLGZU1ZKlUhIwDPwCOMF9/gPgCuC6yLa7vO98GLgPeAy4CTjZbT8VeBJ4CtgD3OK2LwM+AfwEeAB4PzDo9p0F3ABcDOwA3h9TxkHgr4GtwG53zqPdvingRoI3Gd8ITHnf+yZwoct/N/ANYIXbdy9Bo7/HpRPd9rOB24GfAdcA415+CrwDuBPYCXwEEAJxeRzY7/LamXCv3+ry3g3cDbzdbf8lYC/wtFeeI1L8di9yv8EhCfs3Aud4n98GfLdFfqcDP3B5bgVO9e7j+11+ewjeO3YY8Fl37I3ARNl12VI1U+kFsLS4E/A/wJ+7/y91jfzfRrZ90jv+za6BGwL+AngQOMjtey/wmUj+XwI+5hryUeB7XuN+FrAP+FOXXyOmfOuAW4EXOkE5zp3/UCdEZ7rvvtF9Psx975uuoT4WaLjPH3D7JpxgDXnnOR24ywnWEHABsNHbr8BXgOXAGPCwJwJnAde3uc+/CzzPXcMrCDoWv+b2vRK4P+Pv9hZ3Xy4GHnH//563fxfwG97nVcDuhLxe5o5/NcGwyZHAi7z7eJcr+zJgE7AFOMXdpyuBK8qux5aqmWwMziib64CXu/9PBr7tkr/tuvBgVf2Mqu5Q1X2qehGwlEB8DkBEng2sAf5MVX+uqj8laJD98bztqvrPLr+9Mdn8EXCBqm7WgFtUdQeBYNypqv/qvvs54A7gNO+7V6jqFpfvv7PQhRflHcDfq+rtGrhJ/w54qYiMe8d8QFV3quq9BB2DVvktQFW/qqpb3TVcR2BRnpz2+zEcBfwqgTAdAZwLfFpEXuz2H+z2hewCDk4Yh3sbQSfmWlV9WlUfUNU7vP1XuLLvAr4ObFXVGXefvgAc38V1GDXGBM4om28BvyUihwKHq+qdBO6oKbftV/HG30TkL0XkdhfUsJOgV78iIe9xgjG9n4jITnf8xwgsuZD72pTvaAJLLMoRwLbItm0E1kfIg97/vyBo9JMYBz7slfNRAmur0/wWICKvEZHvisijLv81JN+36HfHRGRPmNzmvQTu4Per6pNONP8H+B23fw9wiJfNIcAeVY17AWXSPQ55yPt/b8zn1PfBWFyYwBll8x0CkZomGK9CVR8Dtrtt21X1HgARORn4K4JxuV9W1eUElkFoFUQbz/uAJwjGvpa7dIiqrvSOaffG3/sI3GNRthOIks8YwThfO+LOeR+B63S5lxqqurHD/OYRkaXAF4F/BJ7t7tvXSL5vCzNXvVdVDw6T2/zDNuW4jcCdG3Kc2xZH0j02jK4wgTNKxbnv5oB3ErgmQ6532/zoyRGCMbOHgSERWc9CK+EhYEJEBlzePyFwxV0kIoeIyICIPE9EXpGhiJcDF4rIC1zo+0tE5DACgThWRN4kIkMi8ofAJME4WTseJgjqeK637aPA+SKyEkBElonI76cs40PAUSIynLB/mMCV+zCwT0RewzOWVvj9w0RkWcrzQfC73OvKPCQiJwGvIgiOgWBs7J0icqSIHEEwXvqphLw+AbxVRFa73+hIEXlRhrIYRiwmcEYVuI7AbejP5fq22+YL3DXAfxEEGWwjiB70XYxfcH93iMjN7v+3EDTwmwiCQP4DeE6Gsn2IYPzsGwRRe58gCEbZAbyWoOHeQWBZvlZVH2mXoar+giCQ5gbnkvxNVf0S8EHg8yLyGPAj4DUpy/jfBNbRgyJywPlVdTdwnruOnwFvAr7s7b8D+BxwtyvPESmu4SmCwJg1BFb0x4G3eGNnHyOIeLzVXctX3ba4vL5HEOV5scvrOg60jg0jMxLvEjcMwzCM/sYsOMMwDKOWmMAZhmEYtcQEzjAMw6glJnCGYRhGLen7xZZXrFihExMTZRfDMAzD6JCbbrrpEVU9PO98+17gJiYmmJubK7sYhmEYRoeISHRVoFwwF6VhGIZRS0zgDMMwjFpiAmcYhmEUzuzsLCtXrmTlypXMzs725JwmcIZhGEbhrF+/nk2bNrFp0ybOO++8npzTBM4wDMMonA0bNtBoNAC45557emLFmcAZhmEYPeGYY46h0Wiwd+9e1q9fX/j5TOAMwzCMQgjH3cbHxznttNPYtGkThx9+OI1GgwceeKBwK84EzjAMoyTKCLzoJeG427333svevXtpNBocfPDB7N27l23bthU+FtdTgRORU0Vks4jcJSLvjtn/ThHZJCI/FJFZEbF3QhmGUVvKCLzoJWvXrqXRaDA0FKwpcswxx3DJJZf0bCyuZwInIoPARwhe4jgJvFFEJiOHfR9YpaovIXgxZbNX5TPqy+zsLCeddFLlesh1770b7Vm7di0iUnYxCuOqq65i7969HHvssUxNTXHJJZewevVqrr766vmxuNNOOw1gpJACqGpPEnAicI33+Xzg/BbHHw/c0C7fE044QfuBZrOpIyMj2mw2yy7KomNyclIBbTQaOjMzU3Zx5gnLFabx8fFKlc8onqmpqfm6OTk5Wbvff2ZmRqempmKvK1L/92oRulNEprEngjOAy73PZwKXtjj+UuCChH3nAHPA3NjYWIe3vnfMzMzowMCAAjowMFC7Slx1/AdpcnKy7OKoalAnGo3GAoGroggbxRIKwNjYWG71c2ZmRsfGxrTRaFS6Qx15BvbrYhE44M3Ad4Gl7fLtBwsu7KVZI1YO/oMkIpV46EPRHR4e1tHRURURqx+LmDRehrTC5bc3Ve9Qz8zM6MjIiAbOxP4WuFQuSuAU4HZgNE2+/SBwYS9tfHy8cpbEYqFTK7qI3rAvuGE9iFp0JnKLi7g6EcX3RLSqw81mU4eHh0uvS63ck9HjgD3a5wI3BNwNHAMMA7cAKyPHHA9sBV6QNt9+ELiQKloSi4lmszlvKaV96KPjZHn8ZmEPe2RkZEEZoiJndaSeJI3Ht7Li4lzaSXXYH9crs0MdlmNqaqrtscCc9rPABcPKss4AAAxmSURBVNfAGmCLE7H3uG0bgNe5/2eAh4AfuPTldnn2g8D5PRkbjyuetAPb7UTOF8Q8Ra5VwFFcQ2YiVx9aPf+trDi/UxT1BDWbTW00Gjo6Ojr/eWpqan57WZ2ltBacak0ErojUDwIX7cl0YkkY6fAbiTQ94aT7H22IpqencxOddj3b0C3qn88iLOtBu/H4JCsu2kmOWvp+nr449rJDHdbboaEhXbJkSaY6awLXpwIX9qKiP7ZvSURdVUY64iyhaAMyNjZ2QC8y2kDEuW/C38fv+TabzQV5L126tCOhS9uzjZ7PXJb9T7vx+DRjceFxoXCFaWhoKLae+B3qIl2VUXd+lg68CVyfClwYITQyMrJg+8zMjE5OTs4PBpsll53w3vo907BDET7sYYpaS+3GQ8OHNdogREXHFzo/IGV6enq+Nzs4OLggZRHGODdpp8JqlEP4rPvz3FoJWZwVF9cpigqXL3rR9iauw5b3NSZZlZUegwMOTZGWF1G4NKnKApcmAi/LmJDxDFE3XnjvQgtucnJyQQ83abwr7pi4BsknGqWW5CpqlaINUJZrNbdlf+F7FYaHh+frVZIlH+dmT3JrR/NIGt9Nqut5X2PojWr3DEUpU+AeJ4h+vKdFureIwqVJVRa4NFFEadxlxoEkRYr5D3yz2Zx/qJN+g6hFFg7Qp+l5JgldNHVjwaU5X9bxDiNf2q1SFGeFt6tb0cUJsgRstCpnpx2ksO4l1eFuV2oqU+C+n8cxRaUqC1yWeSA2fSAb4b1tFymW5jeIPviduBGHh4fnhabZbGbqvWYlKbrT3Je9J41lFHaYsgSWtQuW6pRoXW+Xd5L3INp5C68tjTsyjjIF7qCYbb/S7phepSoLXBZs+kDn5OF+iRtb6/Rh7QV+xFqc2JnQ9Ya4wIqkTlboHciyyEBRIpdm5Zy4KSuhoMWJXTftVmkCF/sluLmIwnSSqixwWd0KNn2gc3yByjK+Fc3Dt8T65f77YmdC11tCgRsdHc2loU/KP+/hizSLCvjnjovMbOWyzErVBK40l2Q0VUXg4uaAhBUkiyVgQSedY29sSB6nCxuiJUuW6PLly3V4eLivRLyKRAMp4jqo3Y6dpZ020G3ecfVkyZIlPWuHqiZwf1JEYTpJRQqcL1rteilxroqsA7nhOW1NQqNb0gbAmJXXOXHBSFGrJ5zv1o27O48Ak1Z5txpj61X7UymBq1IqUuCSRCsaFRf1Sfv+7U7cZSZyRl5EXUmhBWfjdt0TJzxxVlE/3NdoPSnCLdqKMoNM2o63lTkm16nAxT34vrUVrahJA6txg6xp5r+1w0SuNUX2ahcDNm7Xfr5jmu8nzWOLWkVVDliKo9fPV5kCtxf4YYt0K302Dy5umZuohRYX7BEVxTzmNrUrp4lcPFlWKjdak+TOrLvQdTPePTPzzHvMWs2xjFumzziQMgVuPEU6qojCpUmdCJxfsaPzOKKp7InXJnIH0m3P24hnsQhd0rhTljmoSa88MjrDxuByErikqKRoQEmVQsXThPQuJsx6K5a6C100ECR6rWmu0yJ288UELieB69eeV9zg9WJ8uMx66x11XBosbtWgTtb6tE5WvpjAJaTh4eHYFdujUWNhwEc/97ziHsR+vI5OKXJOkJFMqykH/SZ2rd6anXSd0euzTlb+lDkGNwH8A/CfwOXAucB4EYXp6ALaRDbG+dn7veeV13vJ+glblLp8QgFIGq9O6mRWxeWftoOUtNZnWP6kVykZnVOmwN0CvA1YDWwDPu7+fgRYWkShsqS0FpxfUeuw1mPc2omdNCJxkaFVFExb4aU6tFsHs10qS+xaWW9RwmtsdX0mcPlRpsD9yPv/++7vELAO+HQRhcqS0o7BhQ15FRvvTmnlOkojWO2mS5Td4w6xNTqrSzQ4KyklCUUvxa4Tyyu8vqjlWodOcpUoU+D+CTjX/X9zZN+WIgqVJVVlLcoyaec6SnIh+ccnrcpStgvUt1StUelf0lh9Ra6Tmce4WR07yVWhTIEbAN4DzAEPAucAb3Yuyu9kOhmcCmwG7gLeHbP/5cDNwD7gjDR5msA9Q1Jvup3LKNqjDSeoRhuiMlxL4WRaFllATZ3J4uLMy2VuUY/VpjSBmz8QngW8DlgPXAz8MbA8w/cHga3Ac4FhgrG9ycgxE8BLgCtN4PIjaQWWdoLVSfRcHsuUFZGXUU2ib+FIWiezW7Hr5+jpxUDpAtf1ieBE4Brv8/nA+QnHfsoErjpkiZ6LHtNN4Ir1uhcnUU9EUp2LdrJaBUxZXao2dRC4M4DLvc9nApcmHNtS4JybdA6YGxsby+cOG23pNnouKobtLErrdRuqC4WrVb1qVe+qME3BSKYogRugD1HVy1R1laquOvzww8suzqJh9erVbNu2jaeeeoprr72WsbExhoaGGBwcnE9LlixhfHyc6elphoeH57eH7N+/PzE99dRTbNu2jVNOOYWhoSHe9a53sXv3bq666qoSr9oom3Xr1vHEE0+wb98+ms3mfL0Skflj9u/fP/+/Xx9DnnrqKR599FFWr17d07Ib5dJLgXsAONr7fJTbZvQhvtjt27dvPj355JP8+Mc/5rLLLptvlKINU1KKNliqysDAABs2bCjxSo0q4YtdtJO1dOlSms3mgvoY1rulS5fyN3/zN2UX3+gxEliHPTiRyBCwhWDC+APAjcCbVPW2mGM/BXxFVf+jXb6rVq3Subm5nEtrlMHs7Cxnn30227dvR1UZGhriwgsvZN26dWUXzTCMAhGRm1R1Vd75DuWdYRKquk9EzgWuIYio/KSq3iYiGwj8r18WkV8HvgT8MnCaiLxPVVf2qoxGuYRWoWEYRh70TOAAVPVrwNci29Z7/99I4Lo0DMMwjK7oyyATwzAMw2iHCZxhGIZRS0zgDMMwjFpiAmcYhmHUEhM4wzAMo5aYwBmGYRi1xATOMAzDqCUmcIZhGEYtMYEzDMMwaokJnGEYhlFLTOAMwzCMWtKztwkUhYjsBjaXXY6cWAbsKrsQOVCX64BqXMsK4JEu86jCdeSFXUv16PY6XqiqI3kVJqSniy0XxOYiXrNQBiJymaqeU3Y5uqUu1wHVuBYRmeu2jlfhOvLCrqV6dHsdIlLIO8/MRVktri67ADlRl+uA+lxLXa4D7FqqSCWvow4uyq57t4ZRZayOG3WnqDpeBwvusrILYBgFY3XcqDuF1PG+t+AMwzAMI446WHCGYRiGcQAmcIZhGEYtqYzAicienPI5VUQ2i8hdIvJub/u3ReQHLm0XkavyOJ9hpKEH9Xu1iNzs6vf1IvL8PM5nGGnpQR3/bVfHfyQinxaRttPcKjMGJyJ7VPXgLvMYBLYArwbuB24E3qiqmyLHfRH4v6p6ZTfnM4y0FF2/RWQLcLqq3i4ifwK8TFXP6rbchpGWIus4cAewDVitqltEZAOwTVU/0Sq/ylhwACJysIjMOpW+VUROd9snROR2Efm4iNwmIt8QkUZMFi8D7lLVu1X1SeDzwOmRcxwC/DZgFpzRUwqu3woc4v5fBmwv+noMI0qBdfww4ElV3eKOuxb4vXblqZTAAY8Dr1fVXwNeBVwkIuL2vQD4iKquBHYSf3FHAvd5n+9323zWArOq+liuJTeM9hRZv/8I+JqI3A+cCXyggPIbRjuKquOPAEMiEs6VOwM4ul1hqiZwAvydiPwQmCG4sGe7ffeo6g/c/zcBEx2e443A57oppGF0SJH1+8+BNap6FHAF8KHui2sYmSmkjmswlvYG4GIR+R6wG9jf7ntVW4vyfwGHAyeo6lMi8mPgILfvCe+4/UBDRI7mmSViPgrcwkJVPwp4IPwgIisITODXF1J6w2hNIfVbRA4HjlPV/+e2/xvwX8VcgmG0pLA2XFW/A5wMICK/AxzbrjBVE7hlwE/djXkVMN7qYFW9D3hp+NlF1bxARI4huClvAN7kfeUM4Cuq+njuJTeM9hRVv38GLBORY90YxauB2wu6BsNoRWFtuIiMqupPRWQp8C7gb9sVphIC5y7qCeCzwNUiciswRxA5kxpV3Sci5wLXAIPAJ1X1Nu+QN2BjE0aP6UX9FpFp4Isi8jSB4J2d4yUYRkt61IavE5HXEgyt/R9V/e+25arCNAEROQ74uKq+rOyyGEbeWP026k5V63jpQSYi8g6CoI8Lyi6LYeSN1W+j7lS5jlfCgjMMwzCMvCndgjMMwzCMIjCBMwzDMGqJCZxhGIZRS0zgDMMwjFpiAmcYhmHUkv8P1gcXYdLGcLsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = ml.read_obs_node()\n",
    "\n",
    "fig, [ax0, ax1] = plt.subplots(2,1, figsize=(6,3), sharex=True, sharey=True)\n",
    "df[ml.obs_nodes[0]][\"theta\"].plot(ax=ax0, marker=\".\", c=\"k\", linestyle=\"\", markersize=3)\n",
    "ax0.set_title(\"Water content at -30 cm\")\n",
    "ax0.set_ylabel(r\"$\\theta$ [-]\")\n",
    "\n",
    "df[ml.obs_nodes[1]][\"theta\"].plot(ax=ax1, marker=\".\", c=\"k\", linestyle=\"\", markersize=3)\n",
    "ax1.set_title(\"Water content at -60 cm\")\n",
    "ax1.set_ylabel(r\"$\\theta$ [-]\")\n",
    "plt.tight_layout()\n",
    "ax1.set_xlabel(\"\")\n",
    "ax1.set_xlim(0,730)\n",
    "ax1.set_xticks([0,365,730])\n",
    "ax1.set_xticklabels([\"Jan-07\", \"Jan-08\", \"Jan-09\"])\n",
    "#plt.savefig(\"../../figures/water_content.eps\", bbox_inches=\"tight\", dpi=300)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[ml.obs_nodes[0]].to_csv(\"../data/wc_30cm.csv\")\n",
    "df[ml.obs_nodes[1]].to_csv(\"../data/wc_60cm.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 9c. Update the soil profile plot with the pressure head"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x10f9d7910>"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANIAAAGoCAYAAAAzV7tRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXiU5dX48e/JZJIhhC2AIouyiKK4oCwigixBUavyaotK1br+qEurbV3esqhYt2pt+9ra2uJSQHHBuqFiQeKCSBFR2RcJIUDCThIgkGWSnN8fMwkBSTJJnsnMPHM+15XLzLPc951hjs8597OMqCrGmMZJiPQAjHEDCyRjHGCBZIwDLJCMcYAFkjEOsEAyxgEWSDFMRI4XkUIR8QRffyYitzZR31eIyJZg/2eJyCoRGRZcN1lEXmmKcUQLC6QoICKDRWShiOwVkTwR+VJE+te1n6puVtVUVS0PoY/JIuIPfvALgv2d24hhPw38Itj/d6raW1U/a0R7Mc0CKcJEpCXwAfBXIA3oBDwMlIShuzdUNRVoDywA3hYROcqYEkNo6wRglcPji1kWSJF3EoCqvqaq5apapKpzVXU5gIgkiMgkEdkkIjtFZLqItAqu6yoiGuIHv4qq+oFpQAegrYjcGDwK/llE9gCTa+pXRJJFpBDwAMtEZENwLNkiMvJo/YnIwOARsEBEllWmgG5igRR53wPlIjJNRC4WkTZHrL8x+DMc6A6kAs82pkMRSQ62uUVVdwcXnwNkAccCj9XUr6qWBI9qAGeqao86+uoEfAg8SuCIey/wloi0b8zfEG0skCJMVfcBgwEFngd2icgsETk2uMm1wJ9UNUtVC4HxwDX1PQoFXSUiBcAWoC9wRbV1W1X1r6papqpFDvZ7HTBbVWeraoWqfgwsAS5pwPijlgVSFFDVNap6o6p2Bk4DOgL/F1zdEdhUbfNNQCKBI0d9zVTV1qp6jKqOUNVvqq3bcsS2TvV7AjAmmNYVBAN5MHBcfQcfzSyQooyqrgWmEggogK0EPoyVjgfKgB1Od33Ea6f63QK8HAzgyp/mqvr7Row16lggRZiI9BKRe0Skc/B1F2AssCi4yWvAr0Wkm4ikAo8TmH0rC/PQnOr3FeAyERklIh4R8YnIsMq/1y0skCJvP4FC/ysROUAggFYC9wTXvwS8DMwHNgLFwC+bYFyO9KuqW4DRwARgF4Ej1H247LMndmOfMY3nqv8rGBMpMRdIInKRiKwTkUwR+W2kx2MMxFhqF7w483vgAiAH+BoYq6qrIzowE/di7Yg0AMgMniQsBV4nUMgaE1ENOTseSZ04/MRhDoEZryoiMg4YB+Dz+fp27Nix6UZnXCcrK2u3qtZ5OVOsBVKdVHUKMAWgR48eumHDFxEekYllIp021b1V7KV2uUCXaq87B5cZE1GxFkhfAz2DZ9uTgGuAWREekzGxldqpapmI/AKYQ+B+mJdUtdaby2bOXNAkYzPxLaamv+sr1mskv7+cnJxCiovrvJPcNJLP56Fz51S8Xs9hy0U6faOq/eraP6aOSPEmJ6eQFi3a0LVrG45yR7hxiKqyZ08+OTn5dOvWqkFtxFqNFFeKi8tp29aCKNxEhLZt2zTqyO/6I1Is10jdu3cnP/9ApIcRNw4cKG7w58VqpCi2Zk0ep5xyUqSHETfWrPmeU05JO2xZqDWSpXamVqmpaXVuc+utt7F69RoAHn/8ycPWDRo0tEF9ZGdnc9ppZ4U4yvrp2vUkdu/eXfeG9WCpXRTr3r07eXmFkR5GnWN46qmnq7Z7/PEnue22O6vWffDBhyH9DUduU1BwkPLyirD8/RUVFeTnHyAhwXfY8sakdq4PpKuuGhzpITTYmjV5pKWl1r1hmKWlpfLZZ58zefKjtGvXlpUrV9G379m88spURIRhwy7g6ad/z7///TZFRUWMGDGc3r1PZcaMaaSmplFYmEdhYSGjR/+Y/PwC/H4/jz46mdGjLz+sj+r27UsBlP/93/tYuPC/dOrUkffee4tmzZqxYcMG7rzzbnbt2k1KSgrPP/93evXqxfvvf8Cjj/6e0tJS2rZNY8aMaRx77LHs2bOHsWOvJzd3K+eeOxARoU2b5j/oc8cO3w8+L1dfHdp75PpAcouH31/H6q37HW3z1I4teOiyk0Pe/rvvlrJq1Xd07NiR884bxpdfLmTw4POq1v/+94/x7LPPsXTp1z/Y1+fz8c47b9KyZUt2797NwIFDuPzyy2qdkVy/PpPXXnuZ559/jquu+ilvvfUO1133U8aNu4N//ONZevbsyVdfLeaOO+7mk0/mMHjweSxa9AUiwgsvvMRTT/2RP/7xKR5++FEGDz6PBx+cyIcfzubFF/9VvzcqBBZIJmQDBvSjc+fAM0v69DmD7OxNhwVSbVSVCRMeYP78BSQkJJCbu5UdO3bQoUOHGvfp1q0rffqcCUDfvmeTnb2JwsJCFi5cxJgxP63arqQk8HTnnJxcrr76WrZt205paSndunUFYP78Bbz99hsA/OhHl9CmzZHP4Gw81weSW2qkX57XKSx9hFq/7NtXREJCYtX2fn8FBQWF5OUVUlZWzt69B6vWHdlmXl4hr776Krm52/n443l4vV769OnD9u15JCWlHnWfgoKDJCZ6q5aXlJRx4MBBdu/eR6tWrfjkk09/0Mftt9/F7bffzsUXX8yCBQt46qmnyMsrpLy8goKCQ+NTVauR6stqpMZLS0ulZctmeL2eqvH4fF5SU32kpaWSmOihVasU0tJS8Xq9tGiRjNfrPWz/8vISOnc+jmOPbcOnn37Gli1baN06paq9o9VIHk9C1fKUlCQqKvx07dqR7t27kZExhzFjfoyqsnz5Cs488wwOHCikV68epKWl8s47b5GYGBjv8OHnM3v2+0yaNJ6PPvoPBQUFjtdINv1tHDVu3C2ccUZfrr32hsOWX3vtWJYs+ZbTTz+b6dNn0KtX6LXZkWbMmMqLL/6LM8/sR+/efXjvvfcBmDx5EmPGjKVv34G0a9e2avuHHprE/Plf0Lt3H95++z2OP/74BvddEzshG8XshGzTaswJWdendm6pkUz42SVCNbAjkqkPu0TImAiz1C6KWWrXtCy1q4GldqY+LLUzJsIskGLKHmC3gz97wjrapUuXMXv2R/Xeb+vWrfzkJ9fUuk1tt1lMm/YyPXueSs+epzJt2sv17r8hrEaKYkfWSGlpTqfhGtYabMGCr1i6dCkDBw4JeZ+ysjJ8vpZMmfJCrWOr6TaL/Px8HnroETIyMhARRowYwZAhw2ndunWdfdslQrVw1yVCxY73UdslSNnZ2Vx00WUMHHgOCxf+l/79+3HTTT/joYceYefOncyYMY0BA/qzePHX3H33PRQXF9OsWTP+9a8pdOvWjaeeeoqioiKWLPma8ePv59JLL+GXv/w1K1euwu/3M3nyJEaPvpypU6fz9tvvUlh4gPLycqZNe4FLL72ClSu/Izs7m+uvv5kDBwK33D/77P8xaNC5P7iEqNKcOR8yatRIevQIPEd01KiRfPXVl4wdW/e1PnYbhQmbzMwNvPnma7z00hT69x/Eq6++wYIFnzJr1vs8/viTvPvuv+nV62S++OITEhMTmTcvgwkTHuStt97gd797kCVLvuHZZ58BYMKEBxgxYhgvvTSFgoICBgw4j5Ej0wH49tulLF++hLS0NLKzs6v6P+aYY/j449n4fD7Wr1/P2LE/Y8mS/9Y43tzcXLp0OfQw3s6dO5ObG/6H8bo+kNyV2jnfR13p0wknnECnTl0pKDjIiSf25JxzziU//wBdunRnw4aN5OUVkpu7jd/+9rdkZWUhIvj9fvLyCiksLKa42F/Vx+zZc3jnnVk8+eQfATh4sIjly9dSWFjM+eefDySRl1d4WNq2b98+7r//flauXInH42HDhg0/2Ka6gwdLKS4uqVpeVFSKakJIKayldrWw1K52taV2+/al0KyZr2qbZs2Sadu2FWlpqezbl4pqBWlpqfzmN39g1Kh07rrrHbKzsxk27ELS0lJJTfXh83mr9vd4Enj33ZmcfPLhF6yuW7eKtLRWVdtVT9v+8pc/c/zxnXj99elUVFTg87UM9n/01K5nz2589tnnVcv37NnJsGFDQ7qK3q7+NhG1d+9eOnUK3C81deqhWbIWLVqwf/+hI8GoURfw17/+ncpzl999tzSkto87rgMJCQm8/PIMystrf/bcqFEXMHfuPPLz88nPz2fu3HmMGnVBQ/6serFAiilOPyjSmfbuv/8exo+fxFlnDaCsrKxq+fDhQ1m9eg19+vTnjTfe5IEHJuD3+znjjL707t2HBx6YXGfbd9zxc6ZNe4Uzz+zH2rXraN68ea3bp6Wl8cADE+jffxD9+w/iwQcnkhaOnPgIrr+y4Yknnoj0MBqse/fudO/ePdLDiBtZWVlkZWUdtuzqq68O6coG1weSXSJkQmWXCBkTYRZIxjjA9dPfbjqPZMLLbqOogdVIpj6sRjImwiy1i2JHpnaPXXIGhXm7HGs/Na09E2cvd6y9I61YsYLt27dzwQX1OyG6bds2xo8fz9SpU2vcZvPmzYwdO5Yvv/zyB+vGjBnDkiVLGDhwIK+99lrI/R6Z2pXW43vHXB9IbrpEyMkgqmwvnA+g3LhxPUuWfMPVV18R8j5lZWWkpfVk1qx/17pdTZcIAUyYcB8HDx7kn/98oV5/35GXCC3PORjyvpbamRplZ2fTq9fp3HjjrZx0Um+uvfYG5s3L4LzzhtGz56ksXhx4WP7ixV9z7rnnc9ZZAxg0aCjr1q2jtLSUBx/8HW+88e+qKxsOHDjAzTePY8CA8zjrrAG8994sAKZOnc7ll1/JiBGjSE+/6LCb9rKzsxkyZARnn30OZ58duJ2jLunpI2jRokWj//6124tC3tb1RyTTOLF2G4WTvrdAOsRNNVI4uO02ikr79hXh95fX6/07skaan5kc8r6uDyQ31Ujh4LbbKCod+dD/UBxZI/35iZUh72s1kmm0aLqNwin7i8vZttcf8vYWSDGkZbtjorK9aLqNAmDIkBGMGfNTMjI+pXPn7syZM7fef9P6nfW7idL1VzbYbRQmVNVvo/h6TyJvbkpm05OX2rdRgNVIJnTVa6TM2bkk5YZ+3s5SO2OOInNnMd3b2axdlVif/t6zZ3+t3/xtnKGqh01/L81qRpfmFSHv7/pAiuXUbuPGvaj6SUtrY8EURqrKnj35tG3bnH79BlPsr+B/v1vGz/p2ZFaIbbg+kGJZ586p5OTks2vX7kgPxfV8Pg+dOwfq0Y27S1CFHu19dex1iAVSFPN6PXTr1irSw4g7G3YFpr57tLcaqUos10gmMuZt8wJJfDv/25D3cX0gxXKNZCJj4evZdCoq5LprBnP92ND2selvY46wYVcx3etRH4EFkjGHUVU27i6p1zkkiIPUzmokUx/7/UJhSQoFW3KZOXNTyPu5PpCsRjL18VVWIaxYz48vPJXzT2pp30ZhTENs3B2Y+u5mqd3hLLUz9fGfXC8e8fLlx1+TUI+LSVwfSJbamfrIeDmLblrCNVefDdgXjRnTINm7S+jaNqne+1kgGRNUUaFsyivhhLb1q48gDlI7q5FMqPaWCsX+FPI25TBzZna99nV9IFmNZEL1VVYhrFzPFRcEpr7BaiRj6m3TnhIAS+2OxlI7E6r/5HpJwMvCeV/zVT3vo4y6QBKRPwCXAaXABuAmVS0IrhsP3AKUA3ep6py62rPUzoTqi9c20tl/kLHBqW+I7dTuY+A0VT0D+B4YDyAipwLXAL2Bi4C/i4gnYqM0rrM5r5Tj0+qf1kEUBpKqzlXVyqcMLgI6B38fDbyuqiWquhHIBAZEYozGnbbkldIlrf7nkCAKU7sj3Ay8Efy9E4HAqpQTXHYYERkHjANo166d1UgmJMXlkHegOQU525g5c3O9949IIInIPKDDUVZNVNX3gttMBMqAGfVpW1WnAFMg8KRVq5FMKNZsK4Jla/nRsJO49Iw2VctDrZEiEkiqOrK29SJyI3ApkK6HnqmcC3Sptlnn4DJjGm1LXmDq2zU1kohcBNwPXK6q1b97cBZwjYgki0g3oCewOBJjNO6zJb8UgC5t3FMjPQskAx8HH4q4SFVvU9VVIjITWE0g5btTVev8jg+rkUwoPt2SRFJCInM/WERDnsXp+m+j2LDhi0gPw8SAcdOz2LSnhDm/PuWw5SKdQvo2iqhL7YyJhJz8Ujo3MK2D6EztHGWpnQnFxp0ptC4rY+bMbQ3a3/WBZNPfpi77i8u5/9vlDO13AlcNPfawdbF8iZAxTSq3IDBj16kRqZ0Fkol7ucGp706trUaqkdVIpi4LdyUCySyZv5T1SQ2bxXZ9IFmNZOqy8T9b8W7dyS0/HUTCEc/gshrJmBBtLSilQyvvD4KoPiyQTNzbWlDKca0aXh9BHKR2ViOZumRubUbX5hWN+qy4PpCsRjK1qahQJi5bxjlndOCqi39we5vVSMaEYveBMvzlSsdGTH1DHByRLLUztdlyIAFoRuaKTGZuWdfgdlwfSJbamdrMXVUA6zbyk0vO5IzOKT9Yb6mdMSHYsc8PQIeW3ka1Y4Fk4tq2vX48CdA2tXHJmetTO6uRTG2+yk6ihcfDW//+slHt2B2yJq5d/2Im+4vLeffOk4+63u6QNSYEO/b5ObaR9RFYIJk4t2Ofn2NaND6QrEYycctfAXuLmrNz01ZmztzUqLZcH0h2HsnUZEteCSxdzYhBJ3JV/7ZH3cbOIxlTh12Fge9qOKZF448nrj8iWWpnarKywAP4+G7RKnYtr2hUW64PJEvtTE1KF+2GrC1ce2X/GmfuLLUzpg679vsRgbTmjT+eWCCZuLWr0E+blES8nobfYl7J1amdqlqNZGr0bWYyvgpx5DPi6kASETImXR/pYZgolZP+Z1IObCdjxpONbstSOxO3ilPa4Tu425G2XH1EUlXSH3050sMwUaioDOYsb87pQ0cy9KphNW43JcRpO1cHkqV2pib7WnWD9D+S89YzZGxdVPcOdbDUzsSloubHAJBycKcj7Vkgmbh0sHkHAJod2O5Ie65O7axGMjV5Z3MSvnzloof+Wet2ViNhNZKp2fpBE/H62pAx6V5H2rPUzsSlopRjaHbAmfoILJBMHFKCgeTQRAO4PLWzGskczX4/zF2RTJ/hozjv6vRat7UaCauRzNEVtOkJw59k88ynKd7+jSNtWmpn4k5RSuAcUrODuxxr09VHJEvtzNF8ut3L8q1w0T2P4fPUvq2ldlhqZ45udZ9xeDsN4suHbnSsTUvtTNwJzNg5l9aBBZKJQ0Up7fE5OPUNLk/trEYyR1KFT5amcPaJHUi/rO7PhtVIWI1kfqg0qQX+S6ex6+PpZGz40LF2LbUzcaUopT0APquRjGm44mAgNSty5hbzSq5O7axGMkf6YmciS3Pgwl89TCiPs7MaCauRzA+tPf0mPN0u4L+Tr6fxT7M7xFI7E1eKUtrhO7jL0SAClx+RLLUzR1q11kdrT+ifC0vtsNTO/NDOS17imG2LyZj2D0fbtdTOxI3yhCRKfa3xFe1xvG0LJBM3ipsFvpXPqaerVufq1M5qJFNd5v4EFqyHQWNv5cQWN4e0j9VIWI1kDpd7/DDodxern72HTQ49z66SpXYmbhQ3awcQlhrJ1UckS+1MdQWbk9hWoFz4u5dC3sdSOyy1M4fbMGginuTWZEy6z/G2LbUzcaO4WTuSw5DWgQWSiSPFzdric/iq70quTu2sRjKVSsphzrLmnDZkJMPHDA15v5ivkUTkHuBpoL2q7hYRAZ4BLgEOAjeq6rd1tGE1kgHgQGpHuPBZct77Oxlb5jveflSmdiLSBbgQ2Fxt8cVAz+DPOOC5CAzNxKiqqxrClNpFZSABfwbuJ/C880qjgekasAhoLSLHRWR0JuZUBlJyUV5Y2o+61E5ERgO5qroskM1V6QRsqfY6J7hs2xH7jyNwxKJdu7ak/81qJAMZ27ys3AYX/+8f8Nbj8BHVNZKIzAM6HGXVRGACgbSuQVR1CjAFoEePHmo1koFDT1ed/+CNYWk/IoGkqiOPtlxETge6AZVHo87AtyIyAMgFulTbvHNwmTF1Km7WNmxpHURZaqeqK4BjKl+LSDbQLzhrNwv4hYi8DpwD7FXVbUdvqao9m/42AKxY4yPNW//PQ1Sndg00m8DUdyaB6e+b6trBpr9Npd2XvETCtsVkTHX2zthKUR1Iqtq12u8K3Bm50ZhYVSGJlPpahzW1i9bpb2McU+JrA4CvOE5qJKdZjWQANhUmMP97GDDmZk5pdUO99nVjjVRvViMZgO0dz4WB97F2ygS27s0OSx+W2hnXK2kWTO3CdAsF1HJEEpErQ9i/WFVnOzgeYxxX4ktDKvx4S/eHrY/aUrvngfeg1qe7nk9gWjoqWY1kAHZlJ1FQ6GFkAz4LTtRIH6lqrc8sEpFX6jOopmY1kgHIHvwQFR4fGa+PD1sfNdZIqnpdXTuHso0xkVbia0NyGKe+IYRZOxHxAD8CulbfXlX/FL5hOcNSOwMwf1kKJ/Y4LqTvjD2Sk9Pf7wPFwAqgot4jiSBL7Ux5QhJF//M6uz6bSca6t8LWTyiB1FlVzwjbCIwJo8qrGpKL88PaTyjnkT4SkQbfH2RMJFUGUlKYAymUI9Ii4B0RSQD8BKbDVVVbhnVkDrAaySzP97B4I5x/y710TKl/ZeJkjfQn4FxgRfAK7JhhNZLZ1P0S6HMrS/94B2tK9oatn1BSuy3AylgLImMASn2tkYpykkr2hbWfUI5IWcBnIvIRUFK5MBamv40p8bUhqWQvQniPA6EE0sbgT1LwJ2ZYjWSyM5PBLw3+HDhWI6nqww0aQRSwGslsHfE0ycX5ZMx4LKz91FkjicjHItK62us2IjInrKMyxiGlya3Cfg4JQkvt2qtqQeULVc0XkWNq2yFaWGoX3yoUPv4uhZPOOZ/0K89tUBtOTn+Xi8jxqroZQEROgDBXbg6x1C6+lSa1pOLSqWybO52MDeG92yeUQJoILBCRzwmcjB1C8JHAxkSzEl+gIkkqDt/5o0qhTDb8R0TOBgYGF/1KVcPzSH9jHFSaHAik5JKCOrZsvNpuNe+gqtsBgoHzQW3bRCOrkeLbd3kelmTD0NvGc4yvYdWIEzXSbODsOvYPZZuIsRopvmWfeCmccTPfPvlzvP4DYe2rtkA6U0Rqu65CgPBed2FMI5Qmt0Iq/CSGOYiglkBSVU/Yew8zS+3i255NSeTta9hDTyrZAyKx1C7ebTp3AuprQ8ak+8Lelz0g0rhWSXIrksJ460R1FkjGtUqbMJBCSu2CTxI6lsOfIrS55j2ig9VI8e2T71I4sdtg0q8Y0OA2HKuRROSXwEPADg49RUiBqH8gitVI8avMk4x/9Gtsz3iNjPXvhr2/UI5IdwMnq2r4nkBujMNKk1sBkFTaNGdoQr3VvGkSTWMcUhlI3kjXSCLym+Cvlbeaf0iM3WpuNVL8Wr3Xw1cb4LybfsPxzRv+XFMnaqQWwf9uDv5Uv9XcbqMwUS3nhBHQ9xcse+bXrD+4M+z91XZlw8MAIjJGVd+svk5ExoR7YMY0RlWNFOanB1UKZbJhPPBmCMuijqV28asoJ4mNu5QLJz+P1PYNX3VodGonIhcDlwCdROQv1Va1BMoaPrSmY6ld/Mrs+wsS25/OJw/8vEn6q+2ItBVYAlwOfFNt+X7g1+EclDGNFbiqoeluTqitRloGLBORVwncMtGLwCTDOlUtbaLxGdMgpUkt8TbROSQIrUa6APgnsIFAQHUTkZ+r6kdhHZkDrEaKX1+vbMbxzStIv7hx//5OP0R/uKpmAohID+BDIOoDyWqk+LX3slfIW5VBxov/apL+QrmyYX9lEAVlEaiTjIlKFQmJlHtToi61WyIis4GZBGqkMcDXInIlgKq+HcbxGVNvpUmBr+5KKm26/9+HEkg+Ald+Dw2+3gU0Ay4jEFhRG0hWI8WnrQcT+Hwt9L3yRk5v07jU3smH6N/UqJFEkNVI8WlP+9NhyMOsnfYoO3evbpI+Q3mI/kkikiEiK4OvzxCRSeEfmjEN408KXCaaVBJdqd3zwH0EpsBR1eXBc0uPhnNgTrDULj4t3JXIsi0w4jdP0MLbuOurnZz+TlHVxXL4BUt2iZCJWpm9xsCpY1n0uxtI0PIm6TOU6e/dwXNHCiAiPwG2hXVUxjSCP6kFif4DTRZEENoR6U5gCtBLRHIJfA3mdWEdlTGN4E9qgbcJp74htFm7LGCkiDQHElQ1Zk7GWo0Un7Izk/GUNfx7Y6tz4jaK39SwHIiNW82tRopP24Y9iddfSMYrjzRZn6Hcan4y0B+YFXx9GbA4nIMypjH8SS1oXti0ZXwot5rPB86uTOlEZDKBi1aNiUpRWSMReMJq9fuPSoPLop7VSPGnXGHOd805ecgFpF81rNHtOXkeaTqwWETeCb7+H2Bqw4bVtKxGij+lSS3h0qls+Wg6GVnh/QLm6kKZtXtMRD4i8CXMADep6nfhHZYxDeNPSgWIytQOVf0W+DbMY3GcpXbxZ1NhAgu+h/5jb6dXq3GNbs++aAxL7eLRrg59YdBEVr3wILn5mXXv4BD7fiTjKn5vZWpX2KT9WiAZV6mqkfxNG0hRmdoFv5PpTqAc+FBV7w8uHw/cElx+l6rOqa0dq5Hiz9ytXtZtV0Y98BwJjXjCaqWYrZFEZDgwGjhTVUtE5Jjg8lOBa4DeQEdgnoicpFrzJb5WI8Wf78+4Bc/xQ/n0gZ81ab/RmNrdDvxeVUsAVLXyqwRGA6+raomqbgQygYZ/p6FxJX9SapPXRxCFRyTgJGCIiDwGFAP3qurXQCdgUbXtcoLLDiMi44BxAO3atSX9b5baxZPszGQSHLryG6I8tROReUCHo6yaSGBMacBAAhfLzhSR7qG2rapTCNw/RY8ePdRSu/iyfegTJJYVk/HKw03ab0QCSVVH1rRORG4H3lZVJXBpUgXQDsgFulTbtHNwmTFVypKa4yva3eT9RmON9C4wHAJPMCLwLYG7CdzGcY2IJItIN6AndjuHOYLfm4rXf6DJ+43GGukl4KXg479KgRuCR6dVIjITWE3g4St31jZjBzb9HW9UYd7SFE4cOMpilysAAA2rSURBVIz0H5/nSJuh1kgS+Iy6U48ePXSkJyfSwzBNpMyTTMbo1+i58mW6f/9O3TuEYMr60m9UtV9d20VjamdMg5R5mwNNf3kQWCAZFzl0eZDVSI6yGim+bCxMYOH30P/aOzmp5e2OtBnV55Gail0iFF92Bm+hWDnlQbYUbGjSvi21M65RFryFItFSO2dZahdfvtyZyIocGHHvH0j1OtOmpXZYahdvNgQfnv/f393YpM/9BkvtjIv4vc3xlBU1eRCBBZJxkTJvCon+gxHp29WpndVI8SU3K5my4gRH/82tRsJqpHizdfBDVHiSyXh1QpP3bamdcY0yb/OITH2DBZJxEauRwsRqpPiyYHkKXbu2J3201UiOshopvhwc/QY7vnyfjCmvNHnfltoZVyhP8KIer6V24WCpXfzY74d5K6D3xWMY9LMrHGvXUjsstYsnB1KPgwv/xoa3n6Noy/wm799SO+MKZYkpABFL7SyQjCuUeSMbSK5O7axGih8rCzwsyYLzfj6RjikVjrVrNRJWI8WT3OOHQb+7+Ob/fs2agzvr3N5pltoZV6hK7cqsRjKmwQ5NNhRFpH9Xp3ZWI8WP4hwv2buUCx6Z6mi7ViNhNVI82dDn5yR0PIeMSTdHpH9L7YwrlCc2i1h9BC4/IllqFz82ZyYjfue+YKySpXZYahdPtg95BFAyZjwYkf4ttTOuUOZtRmJZZGbswALJuER5og9PBAPJ1amd1Ujx48vlKZzQtT3po/s62q7VSFiNFE8OXv4aOxZ+RMaU6RHp31I7E/MqJIGKxGRL7cLFUrv4UFQOHy+DUy64kvOvu9zRti21w1K7eFHcrC1c/DxZ77+IP3teRMZgqZ2JeWWJPgCb/jamMcoTmwGQWFYcsTG4OrWzGik+ZO5PYNF66H/Db+jRwrm7Y8FqJMBqpHix87j+cO54lv/zAbILsiIyBkvtTMwr8wRqpEhOf1sgmZhX7g3USJ6ykoiNwdWpndVI8eHzHYmszoX08c/g8zjbttVIWI0ULzKDX8K84KEbEJydbAiVpXYm5pUn+kgoL4lYEIHLj0iW2sWHfZuT2JmfGJZ/a0vtsNQuXmzq+0sq2vUmY9JtERuDpXYm5pUn+vCUR+6qBrBAMi5QluiL6NQ3uDy1sxopPqxb58MjWI0ULlYjxYc9I57GV7SHjOlPRGwMltqZmFfuSY54ameBZGJeNEw2uDq1sxopPsxflkLX44eSfuW5jrdtNRJWI8WL4tFvsHXB+2SseiViY7DUzsS0CklAPV4Sym36O2wstXO/4uAThHpdcCVDHX6CEFhqB1hqFw9KklvDj14i64OplG2cE7FxWGpnYlq5JxkAT4RTOwskE9PKE6MjkFyd2lmN5H6bDySwcB2c/dNfckqrOxxv32okrEaKB3ntesP5j7DiX4+zdffKiI3DUjsT08o9SYCldmFlqZ37rcj38O1GGHT7QxyXoo63H7OpnYj0Af4B+IAy4A5VXSwiAjwDXAIcBG5U1W/raMtSO5fb2mUo9L+br5+5l+YHtkdsHNGY2j0FPKyqfYAHg68BLgZ6Bn/GAc9FZngmmhxK7UojOo5oDCQFWgZ/bwVsDf4+GpiuAYuA1iJyXCQGaKJHRTCQEiIcSFGX2gG/AuaIyNMEAn1QcHknYEu17XKCy7ZV31lExhE4YtGuXVvS/2Y1kpt9st3L2q0w8oHn8IbhsBDVNZKIzAM6HGXVRCAd+LWqviUiVwEvAiNDbVtVpwBTAHr06KFWI7nb+lOugVOu4vMHr0ciOI6IBJKq1hgYIjIduDv48k3gheDvuUCXapt2Di4zcazCkxR8OGRkRWONtBUYGvx9BLA++Pss4GcSMBDYq6rbjtaAiR/lnqSI10cQnTXS/wOeEZFEoJhgvQPMJjD1nUlg+vumuhqy80jul7cpib37PGH7dw61RhJV509iRYsePXroSE9OpIdhwmh5v19RkNaT8+feGZb2p6wv/UZV+9W1XTSmdsaErNzjxVPuj/QwojK1c4yldu6XnZmMt0wintq5OpDsEiH32zX4YSoSPGS8Mimi47DUzsS0iihJ7SyQTEwr9ySRUGHT32FlNZL7LV3djGN8x5P+I6uRwsZqJPfbd+HfSchbS8ZLf4noOCy1MzGtwuMloaIs0sOwQDKxrSIhiYQomGxwdWpnNZL7fbo0ha7njiR9zPlhad9qJKxGigf+/5lJzhezyFg1I6LjsNTOxCwlAU1ItNQu3Cy1c7fSCpi7FE4aeSXDr7ssLH1Yaoeldm7n96bAZa+Q9Z9XqMj8IKJjsdTOxKyKBC+ATX8b0xhVgWQ1UnhZjeRuu4uF+avh9Ctv4ey2N4alD6uRsBrJ7QpbdIELnmHNzGfJz10Y0bFYamdiVkVC4DiQUBH51M4CycSsykCSivIIj8TlqZ3VSO62sTCBr76HvjfcQ8+WFWHpw2okrEZyuz3tT4chD7PsxUfZvGdNRMdiqZ2JWYdSu8ifR3L1EclSO3dbVeDh2yw45/bJdE6x1C5sLLVzt+2dzoVz7mPJs+NZt29zRMdiqZ2JWRVi09/GNJomeACb/g47q5Hc7avdiazcDEPu/ROtk8LzDHurkbAaye02d78I+oxj4ZN3klyyN6JjsdTOxCwVS+2ahKV27vb5jkTW5sLwSc/h84SnD0vtsNTO7bJOugJOu575D9+MJ8Izd5bamZilVVd/Rz61s0AyMUul8uMbnqsa6sPVqZ3VSO5Wmutl4w5lZBj/ja1Gwmokt9t42vXQ/ZKo+De21M7ELBUPopGvj8ACycQwTUiMmkBydWpnNZK77d2cxJ6CxLD+G1uNhNVIbpdz1m2UdehHxqRbIz0US+1M7ArUSJGf+gaXH5EstXO33dlJHCz0WGoXbpbaudvWfr+iOK0nGZPujPRQLLUzsUslIWpSOwskE7OiKZBcndpZjeRuuRuSkZIEq5HCzWokd9sx8LcUp7Qn49V7Ij0US+1MDIui1M4CycQsFQHC89CT+nJ1amc1krttXJ9McblYjRRuViO5257zHqI8MZmMVyZEeiiW2pnYpSIQJTWSq49Iltq52/rvfSiQfqGldmFlqZ275Q95BFAypj8Y6aFYamdimIBodMzaWSCZmKUkYNPfTcBqJHdbu85HUgKkX2Q1UlhZjeRuBUOfILGsiIxpv4v0UCy1MzEsiq5ssEAyMUuRqJlscHVqZzWSu61a66N5opJ+idVIYWU1krvtG/4UxSV7yZj6WKSHYqmdiW2W2jUBS+3cbeUaHy28SvqP4jS1E5ExwGTgFGCAqi6ptm48cAtQDtylqnOCyy8CngE8wAuq+vsQ+rHUzsX2D/8DpcV5ZEx9ItJDiVhqtxK4EphffaGInApcA/QGLgL+LiIeEfEAfwMuBk4Fxga3NSYqROSIpKprIHDEOMJo4HVVLQE2ikgmMCC4LlNVs4L7vR7cdnXTjNiY2kVbjdQJWFTtdU5wGcCWI5afc7QGRGQcMC74smRK4OhnAtoBuyM9CMesvxuAJXVsVoNQ34sTQmksbIEkIvOADkdZNVFV3wtXv6o6BZgSHMMSVe0Xrr5ijb0fhzj9XoQtkFR1ZAN2ywW6VHvdObiMWpYbE3HRdh5pFnCNiCSLSDegJ7AY+BroKSLdRCSJwITErAiO05jDRGr6+wrgr0B74EMRWaqqo1R1lYjMJDCJUAbcqRr4SjYR+QUwh8D090uquiqErqaE5y+IWfZ+HOLoeyEaJWeGjYll0ZbaGROTLJCMcYBrA0lELhKRdSKSKSK/jfR4mpKIjBGRVSJSISL9jlg3PvierBORUdWWu/L9EpE/iMhaEVkuIu+ISOtq65x7L1TVdT8EJiQ2AN2BJGAZcGqkx9WEf/8pwMnAZ0C/astPDb4XyUC34HvkcfP7BVwIJAZ/fxJ4MhzvhVuPSAMIXlKkqqVA5SVFcUFV16jquqOsqroES1U3ApWXYLn2/VLVuapaFny5iMA5SHD4vXBrIHXih5cUdaph23hS0/sSL+/XzcBHwd8dfS+i7Vo7E6JIXYIVjUJ5L0RkIoFzkzPCMQa3BlJtlxq5gtolWFXqei9E5EbgUiBdgwUSTr8XkS4Gw1RgJgJZBIrIyoKxd6THFYH34TMOn2zozeEFdhaB4tq17xeB+9pWA+2PWO7oe+HKI5KqljXwkiJXaMJLsGLBswSC5ePg/W+LVPU2p98Lu0TIGAe4ddbOmCZlgWSMAyyQjHGABZIxDrBAMsYBFkjGOMACKcaJSFcRqfORYyIyVUQ2ishtDvR5dfAWgw8a25ZbWCDFl/tU9R+NbURV3wBudWA8rmGB5A4eEXk+eDPfXBFpVtcOInJs8Ea3ZcGfQcGj29rg0et7EZkhIiNF5EsRWS8iA+pqN15ZILlDT+BvqtobKAB+HMI+fwE+V9UzgbOBystgTgT+CPQK/vwUGAzcC0xweNyuYYHkDhtVdWnw92+AriHsMwJ4DkBVy1V1b7W2VqhqBYHgytDAdWQrQmw3LlkguUNJtd/LadztMdXbqqj2uqKR7bqaBVL8ygBuBwh+dU6rCI8nplkgxa+7geEisoJAOmjfN9UIdqiOcaqaDZxW7fXTIe63g6M/1KN6WzfW1I85nB2R4sde4BGnTsgCfwfyGz0ql7Ab+4xxgB2RjHGABZIxDrBAMsYBFkjGOOD/A08S5UmKLcRTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 216x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "head = ml.read_nod_inf(times=[730]).iloc[:, 2]\n",
    "ml.profile.loc[:, \"h\"] = head\n",
    "ml.plots.profile()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
